Systems and methods for automated screening and prognosis of cancer from whole-slide biopsy images

ABSTRACT

The invention provides systems and methods for detection, grading, scoring and tele-screening of cancerous lesions. A complete scheme for automated quantitative analysis and assessment of human and animal tissue images of several types of cancers is presented. Various aspects of the invention are directed to the detection, grading, prediction and staging of prostate cancer on serial sections/slides of prostate core images, or biopsy images. Accordingly, the invention includes a variety of sub-systems, which could be used separately or in conjunction to automatically grade cancerous regions. Each system utilizes a different approach with a different feature set. For instance, in the quantitative analysis, textural-based and morphology-based features may be extracted at image- and (or) object-levels from regions of interest. Additionally, the invention provides sub-systems and methods for accurate detection and mapping of disease in whole slide digitized images by extracting new features through integration of one or more of the above-mentioned classification systems. The invention also addresses the modeling, qualitative analysis and assessment of 3-D histopathology images which assist pathologists in visualization, evaluation and diagnosis of diseased tissue. Moreover, the invention includes systems and methods for the development of a tele-screening system in which the proposed computer-aided diagnosis (CAD) systems. In some embodiments, novel methods for image analysis (including edge detection, color mapping characterization and others) are provided for use prior to feature extraction in the proposed CAD systems.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The invention generally relates to diagnostic procedures for thedetection, analysis, assessment, classification, and grading of cancerfrom digitized histopathological images. More particularly, theinvention relates to the use of computer analysis of biopsy samples todetect and classify/grade and score prostate cancer.

2. Description of the Relevant Art

Prostate cancer is the second killer of American men. Current statisticsindicate that 1 in 6 men will suffer from prostate cancer, 1 in 32 menwill die from prostate cancer. Although prostate cancer can be diagnosedat early stages, 80% of men diagnosed with prostate carcinoma are above65 years old.

The risk of developing prostate cancer is commonly attached with thepatient's age. In general, doctors recommend males of age 50 or higherto screen for prostate cancer. However, screening could be recommendedfor males who are at a higher risk of developing prostate cancer even atyounger age, e.g. 40-45. Risk factors generally increase the chances ofdeveloping prostate cancer. The major risk factor is age, since 63% ofprostate cancer cases occur in males of age 65 or older. Furthermore,family history is another major indicator. Men who have a family historyof prostate cancer have higher chances of developing prostate cancer.Considering race, it is reported that African American men have a higherrate of developing prostate cancer compared to Asian and NativeAmericans who have the lowest rate. Considering food diet, it has beenshown that diets that are higher in fats may lead to a higher chance indeveloping prostate cancer.

Biopsies provide the most significant information to urologists byindicating the stage and the grade of prostatic carcinoma using theGleason grading system.

Tissue samples are sampled by urologists from different prostate areascalled “cores” using a hollow needle. When biopsy samples are taken fromprostate tissues, they are prepared for observation under microscope byspecialized pathologists. Cores are sampled from suspicious and nonsuspicious areas in order to have prognostic information on the areaswhere cancerous cells may spread to. Also, sampling different coreswould identify cases where cancer is slowly growing compared to otherareas where cancer is growing more aggressively.

There are different methodologies and protocols that laboratories followin preparing slides of prostate biopsies. In one method, for a singlecore examination, three slides are usually prepared. Each slide may haveserial sections of biopsy slices or multiple slices sections from thesame core. The number of sectioned slices may vary between 2-4 slices.Assuming that four slices are extracted for each slide, the total numberof slices needed to prepare three slides is 12 slices. For serialslices, only the 1st and the 3rd slides are stained using Hematoxylin &Eosin (H&E) or Keratin 903 (K903) prior to pathologists' investigation.Staining is a significant step that clearly differentiates varioustextural components of the biopsy especially basal cells of the prostateglands. Malignant cells do not have these cells and hence, can beobserved clearly. The next step is microscopic investigation where thestained slides 1 and 3 are checked for malignancy. If any suspiciouspatterns are observed, the 2nd slide is then stained using (K903) andobserved. However, if the 3rd slide only shows suspicious patterns,further slices are sectioned from the same core and stained for furtherinvestigation. In case of a negative diagnosis of the current core,slides from other cores of the prostate gland are prepared for furtherdiagnosis.

Pathologists grade cancer tissues by observing biopsy samples under themicroscope to describe their appearance and architectural features byassigning a grade using Gleason grading system. Gleason system ofgrading is a significant histology grading method due to its popularityand usage in many medical institutions as well as in literature. TheGleason grading system was named after Dr. Gleason in 1977, who deviseda system of five grades to describe and differentiate the uniquehistology features of each biopsy image excluding the cytologicalfeatures. The Gleason grading system assigns grades from 1 to 5according five patterns where grade 1 is the most differentiatedpattern, and grade 5 is the least differentiated one. FIG. 1 illustratesthe different architectures of these patterns labeled with theircorresponding grade. Gleason grade 1: cancerous patterns are composed ofgenerally consistent, well separated, and closely packed glands. Thispattern resembles the normal prostate biopsy tissue. Gleason grade 2:tumors of this pattern show less uniformed glands having size and shapevariations, with loose arrangements of the glands and more gaps inbetween. Gleason grade 3: tumor glands in this pattern can berecognizable but with higher variations of shapes (elongated) and mostlysmaller in size. Darker cells' nuclei can be noticed leaving the glandsand invading the surrounding glands. Gleason grade 4: tumors of thispattern can be recognized by the loss of glands formation, and theinability to separate glands from each other. It has an architecturaljump compared to previous grades by having fused and poorly definedstructures. The invading of dark nuclei becomes a dominant pattern inthis grade. Gleason grade 5: tumors of this pattern have no glandulardifferentiation and consist of sheets of single cells. They have a totalloss of architecture with more dark nuclei pattern spread all over theprostate tissue.

Patterns in Gleason grades 1 and 2 are the least common patterns andseldom occur in general. Gleason pattern 3 and 4 are the most commonpatterns, and their existence indicate a worsen prognosis in general.Gleason pattern 5 is less common than the previous pattern and seldomseen as diagnosis usually occurs at an earlier stage of cancerdevelopment.

During a pathologist's examination of biopsy images, they assign twoGleason grades that represent the most two predominant patterns in theimage. The first grade indicates the most dominant pattern, and thesecond grade indicates the second most predominant pattern. The additionof these two grades yields the Gleason score. For instance, a Gleasonscore of 7: 4+3 indicates that prostate cancer tumor of pattern 4 is themost dominant pattern, combined with second predominant pattern of grade3. It is important to note here that a Gleason score of 4+3 is not thesame as 3+4. Although this yields a total score of 7, the first andsecond predominance grades are different at both cases.

During the microscopic observation of Gleason patterns in order toassign the two most predominant grades to the biopsy sample,pathologists may note the existence of a third pattern, or a “tertiarypattern”. Tertiary patterns account for almost 5% of patterns belongingto grades 4 or 5 patterns. Hence, pathologists generally report theGleason score with the tertiary pattern score even if it accounts forless than 5% due to its role in predicting advanced tumor stages at somecases.

Currently, pathologists visually grade prostate biopsy tissues usingGleason scoring system. This approach is very subjective and subject tointra and inter-observer variations. Also, it is a time consuming taskand raises difficulties as far as spatial resolution is consideredespecially in the subgroups of grades 3-5 where further classificationis needed. There is a considerable variation of pathologists' grading ofprostate biopsy images. For instance, at Johns Hopkins Hospital, it wasreported that most “2+2 “biopsies sent for evaluation were actually“3+3” biopsies that had been under graded by a general pathologist.Pathologists experience is another factor that contributes to errors inreporting the correct Gleason grades. Most of these errors result in arepeat needle biopsy procedure for the patient that could have beenavoided if the original biopsy specimens had been reviewed by moreexperienced pathologists. In a study conducted on 3,789 patientscollected from 18 studies, the following were reported: the exactcorrelation of Gleason grade was found in 43% of the total cases, theaccuracy of Gleason grade plus or minus one grade was found in 77% oftotal cases, the under grading of prostate carcinoma using needle biopsyoccurred in 43%, and the over-grading cases accounted for 15% of thetotal cases.

There is a clear need of improved systems and methods for screening,quantitative image analysis, and grading of prostate cancer, providingreliable information on the lesion's extension and localization. Anaccurate diagnosis is a crucial step towards patient treatment selectionand management. Several systems have been proposed for prostate cancerdetection and grading from digitized biopsy slides. The followingpatents and patent applications may be considered relevant to the fieldof the invention:

U.S. Pat. No. 7,761,240 B2 to Saidi et al. discloses systems and methodsare provided for automated diagnosis and grading of tissue images basedon morphometric data extracted from the images by a computer in a2-level procedure (detection and grading). The system proceeds asfollows: (1) Input image, (2) Background removal, (3) Histogram-basedobjects segmentation, (4) Image-level and object-level featureextraction, (5) Feature selection and/or classification. The extractedfeatures include: Fractal dimension of binarized image by severalthresholds, fractal code, variance of symlet4 wavelet coefficients,color channel histograms and tissue structures classification. Themaximum reported accuracy in the detection stage is 96% and 77.6% in thegrading stage. These systems do not have explored the integration amongthe color channels, which could improve the features extraction steps.

U.S. Pat. No. 8,139,831 B2 to Khamene et al. discloses a method forunsupervised classification of histological images obtained from a slidesimultaneously co-stained with NIR fluorescent and H&E stains. Thesystem proceeds as follows: (1) Prostate gland units segmentation, (2)Feature extraction, (3) Principal Component Analysis, (4) Bayesianmulticlass classification (PIN, Gleason grades 1 to 5). (5) Gleasonscore is also computed. The extracted features include: gland boundaryand region description, structural descriptors of glands by usingVoronoi, Delaunay and Minimum spanning tree graph features, as long withtextural descriptors (Haralick features and Gabor filters). In thissystem, only gland centroids are used for graph construction. Some otherimportant structures such as cell nuclei were not explores as a nodes ofof Voronoi, Delaunay and MST graphs. A drawback of this method is thatadditional NIR stain must be used. The most widespread in clinicalpractice is H&E stain.

U.S. Pat. App No. 2010/0098306 A1 to Madabhushi et al. relates tocomputer-aided diagnosis using content-based retrieval ofhistopathological (prostate/breast/renal) image features based onpredetermine criteria and their analysis for malignancy determination.The method classifies a region into normal, cancerous and suspect. Thesystem proceeds as follows: (1) Input image, (2) Feature extraction, (3)Manifold learning (feature dimensionality reduction and classification)for classification and reclassification at higher magnifications Theextracted features include: first and second order statistics of colorchannels, graph features (Voronoi, Delaunay and MST from nuclei centers.The features are standard deviation, average, ratio minimum to maximum,entropy of edges and areas of polygons and triangles, Statistics ofCo-adjacency matrix constructed from the glands centers), texturalfeatures (Haar wavelet feature and Gabor filter response, Sobel,derivatives and Kirsch filter response) and morphology information(nuclear density and gland morphology). The maximum reported accuracy is92.8% in Grade 3 vs. Stroma recognition. However, the classificationbetween Grade 3 vs. 4 is 76.9%.

U.S. Pat. App No. 2011/0243417 A1 to Madabhushi et al. discloses amethod and a system for detecting biologically relevant structures in ahierarchical fashion, beginning at a low-resolution. The relevantstructures are classified using pairwise Markov models (PPMMs). Theinvention provides a fast CAD system for detection/diagnosis of medicalimages. This invention is used for CaP detection in whole mount WMHS,but can be used for other organs. The extracted features include: Graphfeatures (Voronoi diagram, Dealunay triangulation and MST fromlymphocyte centroids) and Gland features and nuclei density. Theperformance of the system is estimated by using AUC of ROC curves. Themaximum reported AUC is 0.812.

U.S. Pat. App No. 2007/0019854 A1 to Gholap et al. discloses a methodand a system that provide automated screening of prostate needle biopsyspecimen in a digital image of prostatectomy specimen. The systemproceeds as follows: (1) Input image, (2) Contrast modification (ifneeded), (3) structure segmentation (lumen, cells, cytoplasm), (4)filtering of isolated cells in stromal area, (5) cell dilation, (6)gland reconstruction, (7) feature extraction, (8) classification usingneural networks. 2-level process (disease detection and furthergrading). The features are: Number of glands, gland size, shape(circularity, elongation), arrangement and destruction, stroma area,lymphocytes presence, among others.

U.S. Pat. App No 2010/0312072 to Breskin et al. discloses a method ofestimating a grade of a prostate cancer from zinc data associated withthe prostate, the zinc data being arranged gridwise in a plurality ofpicture-elements representing a zinc map of the prostate. The methodperforms cancer grading of at least one tissue region based on zinclevels associated with zinc map. In this case, the assessment of cancerseverity is not made from histological quantitative analysis of tissueimages.

It is therefore desirable to have more effective and accurate proceduresfor the detection of prostate cancer using biopsy samples.Unfortunately, all of the above procedures and apparatus sufferdeficiencies in the degree of quantitative analysis and diagnosticaccuracy in the assessment of human or animal tissue images.

SUMMARY OF THE INVENTION

The large number of disparate, differently-scaled Gleason scoringmethods can lead to a variety of misinterpretations of reading actualcancer aggressiveness. We propose generalized scoring methods whichsimultaneously covers all previous scoring methods by relieving thedependency on only two or three of the most dominant patterns, extendingthe range of measured patterns to all possible Gleason patterns, andcombining all observed Gleason pattern grades into one, continuous-scaleGleason score which is bounded by the traditional Gleason Grade limits,1 through 5. These novel scores provide a more accurate, unambiguousmeasure of prostate cancer aggressiveness.

Only few systems have addressed the biopsy whole slide processing. Whilesome of the systems and methods mentioned above are able to classifypredefined regions of interest, they do not apply to a biopsy image as awhole. Certain embodiments disclosed herein map Gleason scores to every,individual pixel within a sample image or group of image slides. Fromeach pixel score, local to global measurements about the entire sampleare made available.

Furthermore, biopsy assessment through integration of both whole slide2-D and 3-D systems is described herein. To the best of our knowledge3-D classification systems have not been investigated for Gleasongrading. We propose a novel approach in cancer diagnosis providing arevised Multidimensional Gleason grading when 2-D and 3-D evaluationresults are provided. Furthermore, systems and methods for cancertele-screening by integrating CAD systems and pathologist experience arepresented in this invention.

To tackle the limitations of tissue cancer imaging systems, the presentdisclosure addresses and overcomes the aforementioned limitations of theprior art by providing methods and systems of automated or partiallyautomated quantitative analysis and assessment of human or animal tissueimages such as liver biopsy, endometrial biopsy, lung biopsy, lymph nodebiopsy, renal biopsy, bladder biopsy, rectal biopsy, skin biopsy, andother serial sections/slides of prostate core images. The methods andsystems are completely general and independent of the process underconsideration.

In an embodiment, a complete scheme/framework for detecting and gradingof cancerous lesions through automated quantitative analysis andassessment of human or animal tissue images is provided. The frameworkdetects, grades, predicts stages of cancers on serial sections/slides ofcore or biopsy images. This output provides information about diseaseseverity, location, and distribution. Specifically, the output of theinvention provides both an overall Gleason grade and a “map” ofcancerous tissue locations. From this, it is possible to determine thesizes and amounts of tissue affected by cancer, the primary samplegrade, secondary sample grades, continuous-scale Gleason score, andhigher-order sample grades. Furthermore, it is possible to determine themost likely distribution of grades within a sample.

The framework of the system consists of a diversity of sub-systems. Eachsub-system could be used separately or in a combined/fused way. Eachsystem use different approaches. For instance, one sub-system is basedon tissue structural and morphology features and other sub-systems arebased on textural characteristics of the tissue.

Also provided are systems and methods for accurate detection and mappingof disease in whole slide images, generating so called “cancer maps” byextracting new features and integrating one or more of theclassification systems mentioned above.

In one embodiment, new methods for image analysis and description ofbiopsy images are based on edge detection, transform based features,color and bit-plans based fractal dimensions, and color mapping aredescribed.

In another embodiment, new methods for image analysis and description ofbiopsy images are based on machine learning, sequential prediction, andinformation analysis.

In one embodiment, a novel “Hard” Gleason score is generated by (seeFIG. 23): Hard Gleason Score=Most Likely Integer Gleason Grade.

In another embodiment, a novel “Soft” Gleason score is generated by (seeFIG. 24): Soft Gleason Score=Most Likely Continuous-Scale Gleason Grade.

In another embodiment, a novel “Hard” Gleason map is generated using themost likely Gleason grade of each sample image pixel (see FIG. 25).

In another embodiment, a novel “Soft” Gleason map is generated using themost likely continuous-scale Gleason grade of each sample image pixel(see FIG. 26).

In one embodiment, a novel “Gleason distribution” is generated using theuncertainty in the most likely Gleason grades of each sample imagepixel.

The scheme/framework based system and sub-systems have many advantages.They are fast, simple, inexpensive, and provides a screening test forcancer that does not give a high percentage of false negative or falsepositive results. The method facilitates early detection of the diseaseand it can be used for screening of large populations. The inventionsystems can be used by pathologists or without pathologists.

In other embodiments, a novel revised Gleason value is generated by (seeFIG. 16): Revised Gleason Value=2D Gleason Grade+3rd Dimensional CoreGrade.

In other aspects, systems and methods for the development of atele-screening system are provided in which the proposed computer-aideddiagnosis (CAD) systems are integrated as biopsy readers according tothe system configuration, pathologists—expert systems level of agreementand diagnosis decision rules.

Embodiments of the computer-aided system can be used for both: (1)assisting pathologists in several stages. For example: localization ofcancer regions in a large digitized tissue biopsy, and grading ofdetected regions, automated quantitative analysis and assessment, (2)Independent assessment without pathologist interaction.

These and other aspects of the present invention will become evidentupon reference to the following detailed description

BRIEF DESCRIPTION OF THE DRAWINGS

Advantages of the present invention will become apparent to thoseskilled in the art with the benefit of the following detaileddescription of embodiments and upon reference to the accompanyingdrawings in which:

FIG. 1 is visual representation of the Gleason grading system;

FIG. 2 is a flow diagram of the Canny's edge detection algorithm;

FIG. 3 is an illustration of the concept for fusion of edge detectedimages;

FIG. 4 is a comparison of the results between one Canny edge detectedimage and the output of the proposed edge detection algorithm;

FIG. 5 is an example of results of tissue structure segmentation usingcolor region mapping in prostate cancer biopsy;

FIG. 6 shows an exemplary flow diagram of 2-D Gleason grading systems;

FIG. 7 shows exemplary morphology- and architectural-based Gleasongrading systems;

FIG. 8 is a block diagram of the first textural-based 2-D Gleasongrading/scoring system.

FIG. 9 is a flow diagram of the second textural-based 2-D Gleasongrading/scoring system using the proposed new Color Fractal Dimension;

FIG. 10 shows two examples of color model transformation for computationof new Color Fractal Dimension;

FIG. 11 is a general flow diagram of the 3-D Gleason grading system;

FIG. 12 presents an example of prostate cancer localization and cancermapping from prostate cancer biopsy;

FIG. 13 illustrates a process to reconstruct 3-D images from 2-D slices;

FIG. 14 shows an example of 3-D reconstruction of segmented lumens andepithelial nuclei from a prostate cancer biopsy image;

FIG. 15 presents a 3-D core grading system;

FIG. 16 shows the steps for computation of “Revised Gleason grade”;”

FIG. 17 is a flow diagram of the integration of CADs with PSA andpatient's information for cancer prognosis;

FIG. 18 shows the structure of the proposed cancer tele-screeningsystem;

FIG. 19 is a method and system for diagnosis of cancer by integratingtwo pathologists and one CAD;

FIG. 20 is a method for diagnosis of cancer by integrating onepathologist and two CADs;

FIG. 21 presents a fully automated method for diagnosis of cancer byintegrating three CADs;

FIG. 22 shows an alternative exemplary classification scheme based onmorphology- and architectural-based Gleason feature selection;

FIG. 23 shows exemplary steps for computation of a “Hard Gleason grade”;

FIG. 24 shows exemplary steps for computation of a “Soft Gleason grade”;

FIG. 25 shows exemplary steps for computation of a “Hard Gleason map”;

FIG. 26 shows exemplary steps for computation of a “Soft Gleason map”;

FIG. 27 shows exemplary results for computation of pattern weights usedfor the “Soft Gleason Score”;

FIG. 28 shows a flow diagram for the implementation of a method forcancer prediction based on histopathology images grading.

While the invention may be susceptible to various modifications andalternative forms, specific embodiments thereof are shown by way ofexample in the drawings and will herein be described in detail. Thedrawings may not be to scale. It should be understood, however, that thedrawings and detailed description thereto are not intended to limit theinvention to the particular form disclosed, but to the contrary, theintention is to cover all modifications, equivalents, and alternativesfalling within the spirit and scope of the present invention as definedby the appended claims.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

It is to be understood the present invention is not limited toparticular devices or methods, which may, of course, vary. It is also tobe understood that the terminology used herein is for the purpose ofdescribing particular embodiments only, and is not intended to belimiting. As used in this specification and the appended claims, thesingular forms “a”, “an”, and “the” include singular and pluralreferents unless the content clearly dictates otherwise. Furthermore,the word “may” is used throughout this application in a permissive sense(i.e., having the potential to, being able to), not in a mandatory sense(i.e., must). The term “include,” and derivations thereof, mean“including, but not limited to.” The term “coupled” means directly orindirectly connected.

Embodiments herein relate to automated quantitative analysis andassessment of human or animal tissue images such as liver biopsy,endometrial biopsy, lung biopsy, lymph node biopsy, renal biopsy,bladder biopsy, rectal biopsy, skin biopsy, and other serialsections/slides of prostate core images. More particularly, but notexclusively, the invention relates to detection, grading, prediction,and staging of prostate cancer on serial sections/slides of prostatecore images, or part of biopsy images, which are illustrated asexamples.

It is to be understood that the present invention is not limited toparticular devices or methods, which may, of course, vary. It is also tobe understood that the terminology used herein is for the purpose ofdescribing particular embodiments only, and is not intended to belimiting. As used in this detailed specification and the appendedclaims, the singular forms “a”, “an”, and “the” include singular andplural referents unless the content clearly dictates otherwise.Furthermore, the word “may” is used throughout this application in apermissive sense (i.e., having the potential to, being able to), not ina mandatory sense (i.e., must). The term “include,” and derivationsthereof, mean “including, but not limited to.” The term “coupled” meansdirectly or indirectly connected. The abbreviation “CaP” refers toprostate cancer.

Described herein is the framework for the development of the cancergrading systems, the functional classification of components and methodsfor image analysis and systems integration.

The Edge Detection Method

Image edge detection is one of the most effective preprocessing toolsthat provide essential image edge information and characteristics. Anedge detector can be defined as a mathematical operator that responds tothe spatial change and discontinuities in a gray-level (luminance) of apixel set in an image. This method can be used in wide areas such asimage segmentation, image categorization, image registration, imagevisualization, and pattern recognition. These applications may vary intheir outputs but they all share the common need of precise edgeinformation in order to carry out the needed tasks successfully. Medicalimages edge detection is an important step towards object recognition ofthe human organs such as brain soft tissues, and other different organs.It represents an essential preprocessing work in medical imagesegmentation. It has also been used for other medical applicationsconcerned with the reconstruction of Computed Tomography (CT) and MRIscans. Conventionally, after image acquisition takes place, various postprocessing algorithms are applied to these images in order to extractmore information. This information can be obtained after edge detection,especially of tissue borders. This process plays a significant rule inrecognition of pathological alternations especially in MRI images.

One of the well-known edge detection algorithms is the Canny edgedetection method. Canny method has proven to be superior over many ofthe available edge detection algorithms and thus was chosen mostly forreal time implementation and testing. It is considered as the modern“standard” in the sense that the validity of all other algorithms isoften checked against it. In Canny algorithm, the Gaussian function isused to smooth the image prior to edge detection. The filtering orsmoothing operation actually services two purposes. The first one is thereduction of the effect of Gaussian noise prior to the detection ofpixel intensity changes. The second purpose is setting the resolution orscale at which intensity changes are to be detected. These factorscontribute effectively towards a precise edge detection method,overcoming the problem of detecting false edges resulted from noisesitting on some pixels.

Canny algorithm has the limitation of being vulnerable to noisedisturbances and impulsive noise. Hence, there are certain limitationsto its application. Also, the traditional Canny operator has a preciseedge positioning but quite sensitive to noise and may detect falseedges, as well as missing some details of real edges in an image. Thisissue is of a great importance in noisy images where many false edgesresulting from noise are detected. Furthermore, the Canny edge detectorcannot detect branching edges to some extent. Because of the existenceof impulsive noise in medical images, traditional Canny algorithmdoesn't perform well since it can't eliminate some noise patterns, thus,generating many false edges. It is difficult to design a general edgedetection algorithm which performs well in many contexts and capturesthe requirements of subsequent processing stages. Consequently, over thehistory of digital image processing, a variety of edge detectors havebeen devised which differ in their mathematical and algorithmicproperties.

Canny operator consists of the main processes shown in FIG. 2. The majorsteps of Canny edge detection algorithm can be summarized as follows:(1) Smooth an image with Gaussian filter to have a noise-free imageprior to applying the derivative function; (2) calculate gradientmagnitude and gradient direction; (3) perform non-maximal suppression toensure the desired edge with one single pixel width; (4) Determine twothreshold values, and then select possible edge points and trace edges.

In one limitation with the Canny method, gradient magnitude anddirection are calculated by pixels within a 2-by-2 neighborhood issensitive to noise and may detect false edges. Another limitation is thefailure to detect branched edges and curves especially in neighborhoodswhere the change in intensity values is inconsiderable. Furthermore,non-maximal suppression can have a negative impact on edge detectionprecision and may miss some edge connected points. However, it becomes avery essential task to detect all of these edges especially in medicalscans where a diagnosis decision is going to be made. The developedalgorithms proposes different modifications to the existing Cannyalgorithm that succeeds in detecting more edge details as well as imagedenoising when compared to the traditional Canny edge algorithm. Also,one might note that the upper and lower threshold levels are fixed atCanny algorithm regardless of the nature of the image. The determinationof such threshold values needs to be an image dependent procedure wherethe characteristics of different images/edge information are utilizedfor an automated selection of threshold levels.

There are a few improved algorithms that were developed based on thetraditional Canny algorithm that has a superior performance to suit ourrecognition purposes. Such proposed algorithms may includereplacements/modifications of different processes in Canny method. Inone embodiment of the present invention, a novel algorithm providesmethods for multiple edge detection. This algorithm has a superiorperformance in localizing and detecting edges that cannot be detectedusing the classical edge detection techniques, including Canny's method.To this end, several modifications to the original Canny's algorithmsare introduced:

-   (1) Improved preprocessing: Gaussian smoothing filter is replaced by    median/adaptive median filter depending on the nature of noise in    the image. This enhances the accuracy of detecting true edges at    areas of interest as well as avoiding the detection of noise edges    that could not be removed using the Gaussian filter.-   (2) Improved calculation of the gradient vector: introducing    different kernel sizes at the traditional 9-pixels neighborhood area    replacing the existing 3×3 kernel. This adds more freedom in    calculating the magnitude at different directions other than 45 and    135 degrees. Also, other algorithms may place certain weights of the    gradient at the direction of interest.-   (3) Improved selection of threshold values: algorithms utilizing    adaptive approach in automating the selection of threshold values    may replace the existing fixed threshold values introduced by Canny.    This improves edge preservation at different areas of the image.-   (4) Improved non-maximal suppression (NMS): Canny NMS process    compares pixels in the n-pixels neighborhood area at eight different    directions in order to determine local maxima. An improved algorithm    considers dividing the area into different sub-areas and applies NMS    at each area before a final decision is made about the local maxima.

In an embodiment, an approach in edge detection is developed based ondynamic range enhancement algorithm. The introduced algorithm succeedsin integrating the concepts of logarithmic ratio of human visual systemwith contrast enhancement in order to develop a new edge detectionmethod. The developed method has the advantage over other traditionalmethods in terms of enhancing, detecting, and segmenting different spotsand shapes having significant dark gray levels. The introduced algorithmis tested on different classes and shows a superior performance comparedto Canny's edge detection operator. Also, it has the advantage of theprecise segmentation and topology/shape preservation of regions ofinterest needed for further measurements and analysis purposes. In theembodiment algorithm, a non-linear image enhancement which is based onHuman Visual System (HVS) dynamic range is integrated with morphologicallog-ratio approach in order to come up with an effective edge detectionoperator that is sensitive to edges of dark areas in the image. Inapplications where edges of interest are represented by dark gray levelssuch as cracks texture, a pre-processing step is needed where spatialfilters such Gaussian/median filters are used to smooth the image priorto applying any sort of enhancement to avoid noise enhancement.

As mentioned previously in an embodiment, a Human Visual System(HVS)-based contrast enhancement and edge detection algorithm may beused to analyze images. The introduced algorithm integrates imageenhancement, edge detection and logarithmic ratio filtering algorithmsto develop an effective edge detection method. Also a parameter ((3) isintroduced to control the level of detected edge details and functionsas a primary threshold parameter. The introduced algorithm functions intracking and segmenting significant dark gray levels in an image.Simulation results have shown the effectiveness of the introducedalgorithm compared to other traditional methods such as Canny'salgorithm in preserving object's topology and shape. The developedalgorithm functions at various applications where measurements andsegmentation of dark gray level spots for classification and trackingpurposes are required such as road cracks, lunar surface images, andremote objects.

The Human Visual System (HVS) is responsible for transferring data intoinformation received by the viewer. The received information can becharacterized by attributes such as brightness, edges, and color shades.Considering brightness, it can be shown that it is not a simple functionof intensity when simultaneous contrasts are considered since it relatesthe perceived brightness of a given region to its background. The reasonbehind this is that our HVS perception is sensitive to luminancecontrast rather than the absolute luminance values themselves. Also, onemight consider eyes sensitivity to white noise added to a uniformbackground more than in an area with a high contrast. Eyes are moresensitive to random noise in smooth areas of an image than in busy areaswith many details. Hence, it is obvious that HVS is not a simplefunction of intensity. Weber's law depends on the ratio of intensityvalues at two points ƒ (x,y) and ƒ_(o) (x,y) more than the difference ofthese values according to the following Weber's law:

$ {\frac{{f - f_{o}}}{f_{o}} = {\frac{\nabla f}{f_{o}} = {\omega \mspace{14mu} ( {Weber}’ s\mspace{14mu} {constant}}}} )$

In another application of parameterized logarithmic image processing(PLIP), images are processed as gray tone function according to thefollowing expression:

g(i,j)=M−ƒ(i,j)

To begin with, the concept of parameterized logarithmic image processing(PLIP) was designed to both maintain the pixel values inside the abounded rang [0,M) where M is the maximum value of the range. Forinstance, M could be 256 for 8 bit gray level images. Also, PLIP modelprocesses images from human visual system (HVS) point of view in a moreaccurate way. The usual arithmetic operations such as addition andscalar multiplication are argued not to be suitable for human visualsystem. The PLIP model considers processing images as gray tonefunctions. The relationship between the traditional gray level functionand the gray tone function is given in as follows:

ƒ(x,y)=M−{tilde over (ƒ)}(x,y)

Where M is the maximum value of the range. The various operationsassociated with PLIP include the PLIP addition, PLIP subtraction, andPLIP multiplications. The operation of interest in the developed edgedetection algorithm is the PLIP addition operation. For two gray tonefunctions a and b, and an arbitrary linear function γ(M)=AM+B forconstants A and B, the PLIP addition is defined as:

${a \oplus b} = {a + b - \frac{a\; b}{\gamma (M)}}$

In order to develop PLIP-based edge detection, Canny algorithmoperations are modified to suit PLIP operations. The gradient magnitudeg and direction θ calculations are replaced by PLIP-based formula asfollows:

Magnitude g=M(1−exp((−√{square root over ( gx ² + gy ²)})/M))

Direction θ=Arc tan( gx, gy)

As an example of this embodiment, gx and gy may be defined in a m×nwindow m,n=1, 2, 3 . . . . The values of m and n may be arbitrary andshould be selected according to the application. For instance, a windowof 3×3 of 9 gray tone pixels may be generated as follows:

$\begin{matrix}f_{1} & f_{2} & f_{3} \\f_{4} & f_{5} & f_{6} \\f_{7} & f_{8} & f_{9}\end{matrix}$ $\begin{matrix}{{g_{x} = {M - {M( \frac{{\overset{\sim}{f}}_{3}{\overset{\sim}{f}}_{6}^{2}{\overset{\sim}{f}}_{9}}{{\overset{\sim}{f}}_{1}{\overset{\sim}{f}}_{4}^{2}{\overset{\sim}{f}}_{7}} )}}},} & {g_{y} = {M - {M( \frac{{\overset{\sim}{f}}_{1}{\overset{\sim}{f}}_{2}^{2}{\overset{\sim}{f}}_{3}}{{\overset{\sim}{f}}_{7}{\overset{\sim}{f}}_{8}^{2}{\overset{\sim}{f}}_{9}} )}}} & {\overset{\_}{g} = {\phi (g)}}\end{matrix}$

where φ(g) is described as isomorphism operator:

${\phi (g)} = {{- M}\mspace{14mu} {\ln ( {1 - \frac{g}{M}} )}}$

This section illustrates the mathematical derivation used in developingthe new edge detector. In current classical edge detection methods, anedge is detected based on a variation of gray level intensity levels.The levels of detected edges details are controlled by the selection ofthreshold levels where a lower threshold value corresponds to more edgedetails and vice versa. In this introduced approach, edges are detectedbased on the difference of adjacent pixels' gray levels instead. Thegray level difference is described for the horizontal and verticaldirections as follows:

Δ_(x) I(i,j)=I(i,j)−I(i−1,j); Δ_(y) I(i,j)=I(i,j)−I(i,j−1)

Such difference in gray levels is applied also using the ratio operationof the same surrounding pixels. The maximum ratio only should beselected when examining rows falling/rising edges, and columnfalling/rising edges.

Rational morphological operations have been investigated in order todesign rational morphological filters:

$\begin{matrix}{{R = ( \frac{{MF}\; 1}{{MF}\; 2} )^{\alpha}},} & {R = \frac{\sum\limits_{i = 1}^{m}{a^{i}( \frac{{MF}\; 1}{{MF}\; 2} )}^{i}}{\sum\limits_{j = 0}^{n}{b^{j}( \frac{{MF}\; 1}{{MF}\; 2} )}^{j}}}\end{matrix}$

Where MF1 and MF2 represent two different classes of filters and a^(i)and b^(i), i=0, 1, . . . m are constants.

Other forms of the developed visual filters deploy logarithmic ratio asthe following examples:

$\begin{matrix}{( \frac{\log \; {MF}_{i}}{\log \; {MF}_{j}} )^{n};} & ( \frac{{\log ( {MF}_{i} )}{{\Theta log}( {MF}_{j} )}}{{\log ( {MF}_{i} )} \oplus {\log( {MF}_{j)} }} )^{n}\end{matrix}$

The basic steps of the new algorithm are the following: (1) Smooth animage with shape dependent filter. (2) Calculate the jointshape-dependent gradient magnitude and direction. (3) Fusion of edgeimages. Some examples of this embodiment are the following:

Step 1: For example the Smoothing operator could be the Gaussian kernel,or

$\begin{bmatrix}{- 1} & {- 1} & {- 1} \\{- 1} & \alpha & {- 1} \\{- 1} & {- 1} & {- 1}\end{bmatrix},\begin{bmatrix}\; & \; & {- 1} & \; & \; \\\; & {- 1} & {- 2} & {- 1} & \; \\{- 1} & {- 2} & \alpha & {- 2} & {- 1} \\\; & {- 1} & {- 2} & {- 1} & \; \\\; & \; & {- 1} & \; & \;\end{bmatrix},\begin{bmatrix}\; & \; & {- 1} & \; & \; \\\; & \; & {- 2} & \; & \; \\{- 1} & {- 2} & \alpha & {- 2} & {- 1} \\\; & \; & {\; {- 2}} & \; & \; \\\; & \; & {- 1} & \; & \;\end{bmatrix},\begin{bmatrix}\; & \; & \; & \; & \; \\\; & \; & {- 1} & \; & \; \\{- 1} & {- 2} & \alpha & {- 2} & {- 1} \\\; & \; & {- 1} & \; & \; \\\; & \; & \; & \; & \;\end{bmatrix}$

Where α=2, 4, or 8.

Step 2: Joint shape-dependent gradient magnitude and direction could bedescribed as:

${{Magnitude}\mspace{14mu} ( {x,y,\theta} )} = {{Max}\mspace{14mu} ( {{\cos \; \theta \frac{\partial G}{\partial x}},{\sin \; \theta \frac{\partial G}{\partial y}}} )}$

Where the gradient magnitude and direction are joined and calculatedusing pixels within an m×n neighborhood, where m,n=1, 2, 3 . . . .

This argument assigns the maximum value of the two functions to themagnitude at a given point, hence, reduces the magnitude value at angleswith low intensity values. Note that the Canny method gradient magnitudeand direction are calculated by pixels within a 2×2 neighborhood issensitive to noise and may detect false edges. Also, a set of definedgradient kernel sets for derivative approximation can be usedeffectively for edge detection, line detection, and consequently forfeature extraction. Such kernels are shape dependent and may vary insizes and dimensions. A kernel size could be of n×m pixels, for exampleof 3×5, 5×7 or 9×5 as shown in the next kernel matrices.

$\begin{bmatrix}1 & \sqrt{2} & 0 & {- \sqrt{2}} & {- 1} \\\sqrt{2} & 2 & 0 & {- 2} & {- \sqrt{2}} \\1 & \sqrt{2} & 0 & {- \sqrt{2}} & {- 1}\end{bmatrix}{\quad{\begin{bmatrix}1 & \sqrt{2} & 2 & 0 & {- 2} & {- \sqrt{2}} & {- 1} \\\sqrt{2} & 2 & {2\sqrt{2}} & 0 & {{- 2}\sqrt{2}} & {- 2} & {- \sqrt{2}} \\2 & {2\sqrt{2}} & 4 & 0 & {- 4} & {{- 2}\sqrt{2}} & {- 2} \\\sqrt{2} & 2 & {2\sqrt{2}} & 0 & {{- 2}\sqrt{2}} & {- 2} & {- \sqrt{2}} \\1 & \sqrt{2} & 2 & 0 & {- 2} & {- \sqrt{2}} & {- 1}\end{bmatrix}{\quad\begin{bmatrix}1 & \sqrt{2} & 0 & {- \sqrt{2}} & {- 1} \\\sqrt{2} & 2 & 0 & {- 2} & {- \sqrt{2}} \\2 & {2\sqrt{2}} & 0 & {{- 2}\sqrt{2}} & {- 2} \\{2\sqrt{2}} & 4 & 0 & {- 4} & {{- 2}\sqrt{2}} \\4 & {4\sqrt{2}} & 0 & {{- 4}\sqrt{2}} & {- 4} \\{2\sqrt{2}} & 4 & 0 & {- 4} & {{- 2}\sqrt{2}} \\2 & {2\sqrt{2}} & 0 & {{- 2}\sqrt{2}} & {- 2} \\\sqrt{2} & 2 & 0 & {- 2} & {- \sqrt{2}} \\1 & \sqrt{2} & 0 & {- \sqrt{2}} & {- 1}\end{bmatrix}}}}$

Step 3: Fusion of edge images. Conventional Canny algorithm uses amethod of non-maximum suppression (NMS) to process the smoothed imageand make sure each edge is of one-pixel width. In this method, a fusionof two edge detection techniques is proposed. The first edge-detectedimage contains the modified kernel set. The second edge-detected imagecontains the same kernel set mentioned previously in addition to usingthe modified gradient magnitude. The fusion of the two edges wouldguarantee the existence of all edges that one image may miss over theother one. FIG. 3 shows the proposed fusion concept.

There are several algorithms that are developed for edges fusion.Considering medical images, the criteria for selecting the fusionoperator of interest must preserve edge information in both images aswell as noise detection avoidance. Hence, the average energy fusionmethod proposed may be utilized. For two edge images G_(L) and G_(H),the fused edge image G_(F) can be defined as:

G _(F,i)=ω_(L) G _(L,i)+ω_(H) G _(H,i)

Where G_(L,i) is the grey-level value at a given pixel location Pi inthe edge image G_(L), G_(H,i) is the grey-level value at a given pixellocation Pi in the edge image G_(H). The average value of thecorresponding coefficients in each edge image ω is defined as:

$\omega = \frac{E_{i}}{E_{L} + E_{H}}$

Where E_(L) and E_(H) represent the low frequency and high frequencyenergy of edge image respectively. Such an average value can judge onthe contents of the fused image, and can reduce the unnecessary detailsthat are contained on one edge image over the other. On the other hand,the energy E_(i) contained in an image at a pixel location P_(j) can bedefined as:

$E_{i} = \frac{\sum\limits_{j = 1}^{N}P_{j}}{N}$

Application to prostate cancer biopsy images: the algorithm is simulatedon a biopsy image having prostate cancer of Gleason grad 3. FIG. 4illustrates simulation results of the algorithm and shows how theproposed edge detection algorithm succeeds in detecting edges moreprecisely compared to canny edge detector.

As mentioned previously, one might consider integrating a properenhancement algorithm as a preprocessing stage for a bettervisualization and hence, an effective detection of different edgedetails.

This combination succeeds in providing useful information about objectsand structural information as well as an enhanced visualization whenvariation of a constant parameter α is considered. In this algorithm,Logarithmic “log” functions replace the morphological filters operation.That is, the deployment of PLIP model while the maximum ratio ofhorizontal adjacent gray level values is investigated. The new rationalexpression becomes:

${{Rx}( {i,j} )} = {\max ( {{\log \frac{( {I( {i,j} )} )^{\alpha}}{( {I( {i,{j - 1}} )} )^{\alpha}}},{\log \frac{( {I( {i,j} )} )^{\alpha}}{( {I( {i,{j + 1}} )} )^{\alpha}}}} )}$

Applying PLIP gray tone function on the row ratio R_(x) yields:

${{gx}( {i,j} )} = {\max ( {{\log \frac{{{int}( {\alpha \; f_{3}} )}{{int}( {\alpha \; f_{6}^{2}} )}{{int}( {\alpha \; f_{9}} )}}{{{int}( {\alpha \; f_{1}} )}{{int}( {\alpha \; f_{4}^{2}} )}{{int}( {\alpha \; f_{7}} )}}},{\log \frac{{{int}( {\alpha \; f_{1}} )}{{int}( {\alpha \; f_{4}^{2}} )}{{int}( {\alpha \; f_{7}} )}}{{{int}( {\alpha \; f_{3}} )}{{int}( {\alpha \; f_{6}^{2}} )}{{int}( {\alpha \; f_{9}} )}}}} )}$

Where ƒ denotes pixels in gray tone; int is the integer operation; a isa predefined constant.

Similarly, the column ratio becomes:

${{gy}( {i,j} )} = {\max ( {{\log \frac{{{int}( {\alpha \; f_{1}} )}{{int}( {\alpha \; f_{2}^{2}} )}{{int}( {\alpha \; f_{3}} )}}{{{int}( {\alpha \; f_{7}} )}{{int}( {\alpha \; f_{8}^{2}} )}{{int}( {\alpha \; f_{9}} )}}},{\log \frac{{{int}( {\alpha \; f_{7}} )}{{int}( {\alpha \; f_{8}^{2}} )}{{int}( {\alpha \; f_{9}} )}}{{{int}( {\alpha \; f_{1}} )}{{int}( {\alpha \; f_{2}^{2}} )}{{int}( {\alpha \; f_{3}} )}}}} )}$

In order to utilize the above expressions to compute the magnitude inLIP edge detection, the isomorphic transformation is applied to have themagnitude components gx and gy. An enhanced magnitude computationconsiders the maximum horizontal and vertical ratios as follows:

${Magnitude} = {M( {1 - {\exp ( {{- \max}\{ {\frac{\overset{\_}{g}\; x}{M},\frac{\overset{\_}{g}\; y}{M}} \}} )}} )}$

Different images were tested in order to verify the functionality of theintroduced algorithm. Images range from medical images to images ofother objects where a different parameter α is selected for each image.The introduced algorithm succeeds in detecting the objects of interestonly by eliminating other unwanted objects that belong to different graytone levels. Also, it has the ability to detect unforeseen edges that noother edge detection algorithm can detect due to the very slight levelof intensity variations. Furthermore, it serves as a trackingapplication to track the development of certain gray level in certainareas.

Most of the traditional edge detection algorithms share one objectivewhich is to detect accurate and more edge details. In an embodiment, theprevious algorithms are utilized in image filtering process. Noisycommunication channels and sensors may contribute to the existence ofimpulsive noise in images. Such images corruption can be eliminated byimplementing image filtering algorithms. Various image filteringtechniques have been introduced but median-based filters have proven toattract much attention due to their ability in preserving edge detailsas well as simplicity in removing impulsive noise.

The following steps summarize the process of image filtering based onedge detection technique: (1) Edge detection processing where edgedetection operator could be any edge detection algorithm or one of thesedescribed previously in this embodiment including shape dependent edgedetection, Noise-resilient edge detection algorithm for brain MRIimages, or logarithmic-ratio based edge detection. (2) Locate edgepixels coordinates in order to avoid such edges during image filtering.(3) Map coordinates to the original/noisy image in order to save graylevel intensity values at these locations. (4) Detect the existence ofimpulsive noise on edges. Detection is based on a rapid change ofintensity values between two adjacent pixels. (5) Apply switching filteralgorithm for image filtering. One filter is a filter of different sizesm×n, e.g. 3×3, 5×5, 7×7 among others square and rectangular kernels,which is applied depending on impulsive noise density. The second filteroperates on noisy edges which is weighted filter with less weight at thecenters of the impulsive noise locations. Switching occurs only wheneverimpulsive noise is detected on edges at locations determined at step (4)previously.

Simulation results of applying the previous algorithm on differentimages were studied. Prior to implementing the developed algorithm,impulsive noise with different densities is added to the image. Thesimulation results show a superior performance of the developedalgorithm in terms of preserving edge details while filtering theimpulsive noise.

In conclusion, an edge preservation filtering technique may be used forimage analysis. The idea behind such a concept is to avoid blurringedges during the filtering of the images, and hence, avoiding missingsome significant branched edges, especially in medical application, e.g.biopsy tissue images. Simulation results showed a superior performancein edge preservation compared to other filtered images where edgesappear to be blurry.

Visual Rational Morphological Filters for Image Enhancement and EdgeDetection

Considering image contrast enhancement, there are various methods forimage enhancement that operate globally on an image such as histogramequalization, gamma correction, logarithmic compression, alpha-rooting,and level or curve-type methods. The drawback of these algorithms isthat they are not always adaptive and in general, they do not take intoconsideration local statistical variations within an image. Also, theydo not match the larger dynamic range the human visual system opposes.Hence, the need to utilize an image enhancement algorithm where pixelsare processed in local neighborhoods is mandatory.An Edge-Preserving Contrast Enhancement (EPCE) is developed based ondynamic range compression. The enhancement part of this algorithm isbased on local mean at each pixel as follows:

${I( {x,y} )} = {\frac{2}{1 + ^{{- 2}{{\tau {({x,y})}}/{\lambda {({x,y})}}}}} - 1}$

Where I(x,y) is the output image, τ(x,y) is either the V component ofthe image in HVS color space or the gray level, and λ(x,y) provides ameasure of the local statistic at a given window for the adjustment ofthe transfer function to the local mean and given by:

${\lambda ( {x,y} )} = {C + {( {M - C} )( \frac{\mu ( {x,y} )}{M( {{Qx} + N} )} )}}$

Where C is an enhancement parameter in the range 0≦C<255, M is themaximum value in that range, Q and N are arbitrary constants. Thisenhancement algorithm functions in an adaptive manner by enhancing darkarea of the image while the light area is preserved. In order to enhanceedges where a significant gray level discontinuity exists, the ratio ofdifference of gray levels is investigated locally to calculate the ratiomagnitude R. An illustrative example is shown for a 3×3 neighborhood asfollows:

$\begin{matrix}{{I( {{i - 1},{j - 1}} )}\mspace{31mu}} & {I( {{i - 1},j} )} & {I( {{i - 1},{j + 1}} )} \\{I( {i,{j - 1}} )} & {I( {i,j} )} & {I( {i,{j + 1}} )} \\{I( {{i + 1},{j - 1}} )} & {I( {{i + 1},j} )} & {I( {{i + 1},{j + 1}} )}\end{matrix}$${{{Rx}( {i,j} )} = {\max ( {\frac{I( {i,j} )}{I( {i,{j - 1}} )},\frac{I( {i,j} )}{I( {i,{j + 1}} )}} )}},{{{Ry}( {i,j} )} = {\max ( {\frac{I( {i,j} )}{I( {{i - 1},j} )},\frac{I( {i,j} )}{I( {{i + 1},j} )}} )}}$${R( {i,j} )} = \sqrt{{{Rx}( {i,j} )}^{2} + {{Ry}( {i,j} )}^{2}}$

The above expression is used to emboss edges details and to yield a newand adaptive dynamic range enhancement as follows:

${Ien}_{({i,j})} = {\frac{\alpha ( {R( {i,j} )} )}{1 + ^{{- 2}{{\tau {({i,j})}}/{\lambda {({i,j})}}}}} - 1}$

Where len_((i,j)) is the enhanced pixel, R(i,j) is the ratiocoefficient, a is an arbitrary function operating according to the valueof R. The proposed formula succeeds in enhancing and highlighting edgeswhere gray levels transitions are significant.

Recall from the log-ratio approach for morphological filters. The effectof such ratio is to give more weight to edges of interest that need tobe enhanced to the highest dynamic range of the HVS based on the localstatistics as follows:

${I_{ed}( {i,j} )} = {{\beta ( {\log \frac{\tau ( {i,j} )}{F( {i,j} )}} )}I_{{en}{({i,j})}}}$

Where β is a normalization coefficient controlling the enhancementsaturation level and acting as a primary edge detection threshold,F(i,j) is local a filtering operation based on local neighborhoodstatistics. A further threshold stage is introduced as an adaptiveprocess to reflect local statistics by deploying standard deviation inlocal windows. Several thresholds may be used. Particularly, we use athreshold T(i,j) for illustrative purpose which is calculated as follows[16]:

T(i,j)=α·tm(i,j)+C·σ(i,2,5)

Where α·tm(i,j) is the local alpha-trimmed mean, C is a global constant,and σ(i,j) is the standard deviation. Once potential edge gray levelsare detected, binarization process takes place. A summary of thedeveloped algorithm is provided in the table below:

Process Algorithm Pre-processing Gaussian/Median/Order-statisticsfiltering Ratio magnitude calculation R(i, j) = {square root over(Rx(i, j)² + Ry(i, j)²)}{square root over (Rx(i, j)² + Ry(i, j)²)}Adaptive enhancement${I_{en}( {i,j} )} = {\frac{\alpha ( {R( {i,j} )} )}{1 + e^{{- 2}{{\tau {({i,j})}}/{\lambda {({i,j})}}}}} - 1}$HVS-based edge detection${I_{ed}( {i,j} )} = {{\beta ( {\log \; \frac{\tau ( {i,j} )}{F( {i,j} )}} )}{I_{en}( {i,j} )}}$Threshold T Adaptive thresolding Example: T(i, j) = α · tm(i, j) + C ·σ(i, j) Binarization${I_{ed}( {i,j} )} = \{ \begin{matrix} 1arrow{{I_{ed}( {i,j} )} \geq T}  \\ 0arrow{elsewhere} \end{matrix} $

The architecture and characteristics of dark nuclei is an essentialfeature in deciding the grade of cancer besides other glandularfeatures. In this case, the proposed algorithm operates on the image andsucceeds to segment dark nuclei. This result can be used as an input tofeed classification systems of interest such as fractal analysis basedclassifier for features extraction.

In this embodiment, we also present an algorithm for visualizing edgesin colored images as well as edge detection using the previouslydeveloped logarithmic ratio approach. The developed visualizationalgorithm has a superior performance in highlighting more unforeseencolor edges than the standard color-grayscale conversion method canpresent. Also, by integrating the proposed algorithm with alogarithmic-ratio based edge detector operator, the developed algorithmoutperforms the standard edge detection operators in gradient coloredimages where color boundaries transitions are hard to detect.

Most of edge detection operators operate on grayscale images since colorimage components impose some sort of complexities in detecting edges. Ingrayscale images, an edge can be detected by investigating thediscontinuity in gray levels. However, in color images there are noclear and defined color edges although color images carry moreinformation and therefore have the potential of having more hidden andaccurate edges. The limitation factor in detecting all the edges in graylevels images is due to the fact that different colors have differenthue values but they share the same intensity value. Hence, an effectiveapproach to convert color components from several color models such asRGB components into gray level scale is needed.

Considering color space, there are different representations of threedimensional color spaces having different characteristics. One set ofcolor spaces uses Cartesian coordinates for point's representation inspace. Examples of this type of sets are the three primary illuminationcolors Red, Green, and Blue (RGB), the complementary colors Cyan,Magenta, and yellow (CMY), and the opponent color representationluminance, in-phase and quadrature (YIQ). Another space representationutilizes polar/cylindrical coordinates such as Hue, Saturation, andLightness (HSL) and Hue, Saturation and Value (HSV). Generally, eachcolor space can be converted to another color space representation.

For color images edge detection, previous researches in this areautilized vector order statistics, distribution of pixel colors andfusion techniques of edges detected at different color spaces. Hence,more focus was given to analyzing and extracting color informationrather than investigating new edge detection operators.

Considering real world images, one might consider the precision of edgemaps as a result of human visual system (HVS) perception. One of thewell-known image processing frameworks that incorporate thecharacteristics of HVS in image enhancement as well as edge detection isthe Parameterized Logarithmic Image processing (PLIP) model. In PLIPmodel, a new set of arithmetic operations are introduced taking intoconsideration the characteristics of human visual system.

Furthermore, the concept of utilizing ratio operations in morphologicalfilters design that oppose some characteristics of HVS was introduced.In previous embodiment, the logarithmic image processing framework andthe concept of rational filters design were integrated to produce a newedge detection algorithm based on HVS model. The proposed methodprovided a new approach in edge detection focusing on the class ofimages where the variations and discontinuities in gray levels follow agradient pattern. In this embodiment, we introduce an algorithm forvisualizing unforeseen edges represented by any color model. An exampleis presented taking the representation of an image in RGB color space.

RGB color model consists of three independent color planes representingthe three chromaticity components: Red, Green, and Blue. In order tospecify a pixel RGB value, each color component is specifiedindependently within the range 0 to 255. RGB color model is used in manyapplications such as computers, television, and electronic systems.Conversion of RGB color space to grayscale is the common preprocessingstep prior to applying edge detection. The concept of converting RGB tograyscale depends on assigning weight for each color component based onHVS and its perception of colors. The common weighting vector used is[L, M, N], where L+M+N=1.

Some of the common edge detection techniques for color images utilizefusion technique of edges detected at different color space usingtraditional edge detectors. The first fusion method finds the gradientand the edge map at each of the three RGB levels, then fuses theresultant edge maps output. The second method fuses the gradients ateach RGB component, and then finds the edge map of the final fusedgradient.

The developed algorithm focuses on extracting edge details that cannotbe highlighted for detection using the standard assigned HVS vectorweight. To illustrate the operation of the algorithm, the RGB colormodel is considered. Given a tri-chromatic color image, one mightconsider designing a rational morphological filter that takes intoconsideration the assigned HVS weight for each component intri-chromatic space. For instance, for a RGB image let and Fn and Fmrepresent two different filters as follows:

Fn=(L*R,M*G,N*B)

Fm=[βR log^(αR)(L*R),βG log^(αG)(M*G),βB log^(αB)(N*B)];

[βR log^(αR)(w ₁ R−w ₂ G),βR log^(αR)(w ₁ R−w ₂ B)];

[βG log^(αG)(w ₂ G−w ₁ R),βG log^(αG)(w ₂ G−w ₃ B)];

[βB log^(αB)(w ₂ B−w ₁ R),βB log^(αB)(w ₂ B−w ₃ G)];

Where β and α are constants assigned to each color component, w_(i), w₂,w₃ represent different assigned weight for R,G,B color componentsrespectively. The ratio of different morphological operations may takeany of the following forms:

${\beta \; 1( \frac{\max ({Fn})}{\max ({Fm})} )};{{\beta 2}( \frac{\min ({Fn})}{\max ({Fm})} )};{{\beta 3}( \frac{\min ({Fn})}{\min ({Fm})} )};{{\beta 4}( \frac{\max ({Fn})}{\min ({Fm})} )};$

where β is normalization constant

The above expressions yield successful results in investigating theintensity of each RGB component at each pixel taking into considerationthe HVS perception characteristics. The result of carrying out thesedifferent operations produces a gray level image having more edgedetails than the standard translation of RGB to grayscale.

Image processing algorithms presented herein highlight and visualizedifferent information that can be obtained from variety of imagesincluding biopsy images as well as real world images. The developedpreprocessing algorithms succeed in detecting, segmenting andvisualizing different image details compared to the state of the artcurrent image processing algorithms.

Method for Carcinomas Color Region Mapping

Referring now to another embodiment of the invention, a method forfinding carcinomas dominant colors is provided. The goal of this methodis to reduce the number of distinct colors or define the most populatedcolor in the biopsy image. The resulting quantized image should be asvisually similar as possible to the original image for the furthertissue structure objects segmentation. For example we may map a biopsyimage containing millions of colors into an image containing a fewcolors, which are closely tied to tissue structures in H&E stainedneedle biopsy images are also included in the color map. Additionalcolors of interest can be added manually to keep information for furtherstudies.

Before performing color region mapping, color standardization in digitalhistopathology images is performed. Various algorithms me be used tomake the appearance of the image similar to those in the trainingdatabase and tackle the problem of color constancy. Examples of colornormalization algorithm include: Reinhard's method and its modifications[18], color deconvolution, histogram matching, among others.

The general algorithm for constructing color region mapping comprisesthe following steps: (1) Survey the colors which compose biopsy images.(2) Determine the ‘best’ set of colors for a biopsy image by usingalgorithms for localization of local maxima in 2-D color histograms. (3)Generate color table. (4) Assign each color to a class for segmentationpurposes. (5) Map images into the reduced color space.

For instance, to construct the color region mapping shown in thefollowing table, instead of adopting an absolute set of color featuresfor our segmentation algorithm, the 2-D projections of 3-D RGB colormodel of a large variety of prostate needle biopsy images of differentgrades are computed to construct 2-D histograms and generate a goodapproximation of color region distribution in those images. This is animportant step since the variations in colors of stained imagesconstitute an aspect to take into account before applying segmentationprocedures based on their color composition, since these variations arealways present even if there are defined laboratory controls andprotocols to prepare tissue core samples. Based on the calculated 2-Dhistograms, a color region mapping is proposed taking the colors thatappear more frequently in our image set and that represent any object ofinterest in the evaluation of histopatholy images.

In more detail, still referring to the carcinomas color region mappingmethod, the most prominent colors may found by using an automaticalgorithm for local maxima location in 3-D surfaces.

The procedure explained above is valid for finding representative colorsof any image dataset. In this embodiment, the resulting region mappingfor prostate cancerous class of test images is presented in the nexttable.

COLOR COMPONENTS STRUCTURE i COLOR R G B (class) 1

255 255 255 Lumen and non- tissue areas 2

205 146 178 Epithelial cytoplasm 3

219 133 168 Stroma 4

225  77 109 5

159 167 203 Blue mucin 6

 90  37  90 Epithelial nuclei

It is important to note that any color model can be used to apply theabove method. In this case, color model transformation must be performedbefore constructing 2-D projections.

This embodiment also provides a method for tissue structuresegmentation. Segmentation procedure is carried out through pixelclassification. This automatic process consists of minimizing thedistance between each pixel and the colors considered in the color map,and assigning to each pixel the class corresponding to the nearestcolor. To classify a pixel, an appropriate distance metric may be used.As an example Euclidean distance is:

d _(i)(m,n)=√{square root over ((r _((m,n)) −r _(i))²+(g _((m,n)) −g_(i))²+(b _((m,n)) −b _(i))²)}{square root over ((r _((m,n)) −r_(i))²+(g _((m,n)) −g _(i))²+(b _((m,n)) −b _(i))²)}{square root over((r _((m,n)) −r _(i))²+(g _((m,n)) −g _(i))²+(b _((m,n)) −b _(i))²)}

In the equation above d_(i)(m,n) is the distance in RGB color model frompixel (m,n) to the i'th color region defined in the color map table.FIG. 5 shows the structure segmentation results based on colorinformation of an image that belongs to prostate cancer, Gleason grade 3pattern.

Method for Pattern Detection

The novel Gleason grading concept exploits the ability of an optimalpredictor/compressor to compress (in an informational sense) itspredicting class optimally. In this embodiment, an implementer selectsan appropriate machine learning technique and an appropriate predictiontechnique (which may be one in the same) to train a set of predictors onimages or pre-processed images containing a plurality specific Gleasonpattern grade. By pre-processed we imply that operations may have beenperformed on the image (such as architecturalsegmentation/classification, textural segmentation/classification,feature segmentation/classification, denoising, scaling, etc.) prior topredictor training. Subsequent use of the word “image” with regard topredictors refers to an image that has already been appropriatelypre-processed.

Each predictor should be trained on a specific type of Gleason pattern.Possible techniques include: maximum likelihood estimation, Bayesianestimation, Kalman filtering, prediction by partial matching, contexttree weighting, and others. Paired with an appropriate entropy encoder,each predictor/encoder pair forms the best compressor for its respectivepattern type compared to the other predictor/encoder pairs. Hereafter,predictor/encoder pairs shall be referred to as pattern compressors(PCs) for simplification. Scoring of a biopsy sample image (hereafterreferred to as the query or query image) amounts to compressing eachquery pixel with each PC. Smaller compressed sizes per pixel imply moresignificant membership to the pattern class of a PC. Selection,combination, or other analysis of the all PC compression rates can thenbe used to estimate the grade of a particular pixel or the query imageon a whole.

As an exemplifying embodiment, universal detector/classifiers (UDCs)constructed from a variable order Markov model which utilizes aprediction by partial matching (PPM) algorithm for prediction and anarithmetic encoder for compression can be used as PCs. Each UDC trainson a given Gleason pattern type such that each is the best UDC datacompressor of that pattern type on average. Grading of a query sampleinvolves selection, combination, or other analysis of the classifier'srelative compression performance. Optionally, a PC optimized to thequery image can be used for comparison purposes.

Method for Hard Gleason Scoring & Mapping

A “Hard Score” may be assigned to a biopsy image to directly correspondto the single most-likely pattern observed within a query image. This isanalogous to a pathologist's choice of grade, where a hard decision mustbe made between possible Gleason grades. This is achieved by choosingthe predictor/encoder pair which best compresses a query image onaverage. The hard score of the sample is the grade of the pattern forwhich the predictor/encoder pair was optimized (see FIG. 23).

The same hard grading method may be used on a pixel by pixel basis,assigning a hard Gleason grade to each pixel individually depending onwhich trained UDC compresses it best. This results in a map of pixelgrades which can distinguish between regions of different Gleasonpatterns within a single biopsy image (see FIG. 24). It should be notedthat the average hard map grade is not necessarily equal to the hardgrade as determined above. For instance, if patterns 3 and 5 are bothprevalent within a query image, then the average hard score is likely tobe near 4 despite grade 4 Gleason patterns being absent within theimage. The hard score, however, will be 3 or 4 depending on whichpattern is most prevalent.

Method for Soft Gleason Scoring & Mapping

A biopsy image may also be assigned a “Soft Score,” here denoted s. Asoft score can better characterize a biopsy image which is a member tosome extent of one or more predictor/encoder pair classes. One way togenerate such a score is to assign membership values, or weights w_(g),to each possible class found within an image and to compute the weightedmean. Thus, the score is a weighted combination of all Gleason gradeswhere the grades G_(i) range from 1 to 5 and the weights w_(i) aredetermined from a query image:

s=w ₁ ×G ₁ ⊕w ₂ ×G ₂ ⊕w ₃ ×G ₃ ×w ₄ ×G ₄ ⊕w ₅ ×G ₅

In the above equation, the x symbol indicates a scaling operationdefining how much of each pattern grade is within a given sample. ⊕ is ageneralized combination operator defining how the weighted grades shouldby joined into one, global Gleason score. Other embodiments may include,but are not constrained to:

$\mspace{20mu} {s = \sqrt[{w_{1} + w_{2} + w_{3} + w_{4} + w_{5}}]{G_{1}^{w_{1}}G_{2}^{w_{2}}G_{3}^{w_{3}}G_{4}^{w_{4}}G_{5}^{w_{5}}}}$$\mspace{20mu} {s = \frac{w_{1} + w_{2} + w_{3} + w_{4} + w_{5}}{\frac{w_{1}}{G_{1}} + \frac{w_{2}}{G_{2}} + \frac{w_{3}}{G_{3}} + \frac{w_{4}}{G_{4}} + \frac{w_{5}}{G_{5}}}}$  s = arg  max (w_(g), G_(g)), g = 1, 2, 3, 4, 5  s = arg  min (w_(g), G_(g)), g = 1, 2, 3, 4, 5s = PLIP(w₁ × G₁) + PLIP(w₂ × G₂) + PLIP(w₃ × G₃) + PLIP(w₄ × G₄) + PLIP(w₅ × G₅)

One embodiment is the form the above equation takes by using standardmultiplication as the scaling operation and addition as the combinationoperation. The equation then becomes linear combination of patterngrades:

$s = \frac{{w_{1} \times G_{1}} + {w_{2} \times G_{2}} + {w_{3} \times G_{3} \times w_{4} \times G_{4} \times w_{5} \times G_{5}}}{w_{1} + w_{2} + w_{3} + w_{4} + w_{5}}$

While traditional binary and tertiary scores prove a score that is thesum of the weighted grades, the above equation normalizes the score bythe sum of the weights to ensure that the final measure is a scoreranging continuously between 1 and 5.

The relation to the traditional binary Gleason score can be described asfollows: the traditional binary score assigns a weight of 1 to the twomost prevalent pattern grades and a weight of 0 to the three remainingpattern grades without normalizing the sum by the sum of the weights.The modified and previously-mentioned continuous-scale scores can besimilarly defined through appropriate definition of the pattern weightsw_(i) and the removal of the normalization factor. The benefit of anormalized system is that it provides a continuous, threat level graderanging from 1 to 5 in correspondence to the traditional Gleason gradingsystem.

Often, only grades 3, 4, and 5 patterns are considered in medicalpractice. Thus, if we have available to us 3 predictor/encoder pairs,each optimal for a grade 3, 4, or 5 Gleason pattern, thenpredictor/encoder pair's respective pattern grade is scaled by itsrespective weight for averaging:

$s = \frac{w_{3} \times 3 \times w_{4} \times 4 \times w_{5} \times 5}{w_{3} + w_{4} + w_{5}}$

The weights w₃, w₄, and w₅ are the membership values determined throughstatistical analysis of the compression rates Y_(PC) obtained for eachtrained PC on a query image combined with compression rates obtained bya fourth PC, PC₀, trained on the query image itself. Compression by PC₀generates a more optimal self-compression rate which is the amount ofinformation the query image provides about itself.

For one embodiment utilizing analysis of variance (ANOVA) as astatistical comparison tool, confidence interval estimates Ŷ_(PC) forthe four PC mean compression rates are t-distributed around theirrespective measured mean compression rates, Y _(PC), parameterized bythe degrees freedom in the sum squared error (which depends on bothsample size and the number factor levels in question). When the samplesize is large enough, the t-distributions are practically

${\sigma^{2}\{ {\hat{Y}}_{UDC} \}} = \frac{MSE}{n}$

normal, centered at Y _(PC) with variance where n is the number ofsamples (which is just number of query image pixels) and MSE is the meansquare error as determined through 4-level, single factor ANOVA testingof the four PCs' compression rates.

FIG. 27 gives an example of the four PC estimate intervals asoverlapping normal distributions centered at their respective meancompression rates. The total area of the intersection between the PC₀distribution, Ŷ₄, and a particular Gleason pattern-trained PCdistribution, Ŷ_(g), is a measure of how similar the curves are to oneanother and therefore provides a membership value of the query image tothat particular PC and can be used as a weight. For the example in FIG.27, the amount of overlap between the compression rate estimates for PCs3, 4, and 5 with that of PC₀ are 2.55e-10, 3.51e-4, and 1.23e-4,respectively. Then, the soft score s is computed:

$s = {\frac{{2.55e} - {10 \times 3} + {3.5\; 1\; e} - {4 \times 4} + {1.23e} - {4 \times 5}}{{2.55e} - 10 + {3.5\; 1\; e} - 4 + {1.23e} - 4} = 4.26}$

A 3-level, single factor ANOVA test provides information for whether ornot the three Gleason pattern-trained PCs compress the query imagesignificantly differently. If they do not, then the soft score isunreliable and the image either belongs to all classes equally or is notlikely a member of any of the classes (at least as defined by thetraining set). In this case, the soft grade returns “CANNOT BECLASSIFIED” which may prompt a pathologist to manually inspect an imageand possibly add it to a training set for the proper PC.

Like hard Gleason mapping, soft Gleason mapping assigns each pixel anGleason grade dependent on the compressibility of a query image pixelI(i, j) by each PC. However, rather than selecting the best-compressingPC class, soft grading factors in the performance of all PC classes togenerate an average grade. Thus, the soft grade of each pixel is aweighted mean of each Gleason class. If only grades 3, 4 and 5 areconsidered:

${s_{ij} = \frac{{w_{3,{ij}} \times 3} + {w_{4,{ij}} \times 4} + {w_{5,{ij}} \times 5}}{w_{3,{ij}} + w_{4,{ij}} + w_{5,{ij}}}},$

or if a PC of each pattern class is available:

$s_{ij} = {\frac{{w_{1,{ij}} \times 1} + {w_{2,{ij}} \times 2} + {w_{3,{ij}} \times 3} + {w_{4,{ij}} \times 4} + {w_{5,{ij}} \times 5}}{w_{1,{ij}} + w_{2,{ij}} + w_{3,{ij}} + w_{4,{ij}} + w_{5,{ij}}}.}$

Computation of the weights w_(g,ij) may differ substantially from themethod used to compute the global weights w_(g) because global imagecompression statistics may neither be available for nor relevant to thelocal grade. However, PC predictions are based on a pixel's localcontext and are locally applicable. In fact, a PC estimates the locallikelihood function defining the probability of observing a particularpixel's value under the assumption that local the pixels are members ofthe PC's class. The observed likelihood models produced by each PC afterobservation of a pixel can be compared to the optimal likelihoodfunction as measured by the a query image's self-PC, PC₀, and measuredas an error.

A convenient measure of error (and the error for which is used to proveuniversality) is the Kullback-Leibler divergence. The KL divergencemeasures the amount of extra bits produced using one probabilistic modelinstead of a more accurate probabilistic model. In this case, we comparethe likelihood function for a subsequent observation generated by a PCconditioned on a specific Gleason pattern to the true likelihoodfunction generated by the PC conditioned on the query image I afterobservation of pixel I(i,j). The KL divergence of PC_(g) to PC₀ afterpixel I(i, j) is given by

$ɛ_{ij}^{{PC}_{g}} = {\sum\limits_{a}^{\;}\; {{L^{{PC}_{0}}( {a{{I( {1,1} )}\mspace{14mu} \ldots \mspace{14mu} {I( {i,j} )}}} )} \times {\log_{2}( \frac{L^{{PC}_{0}}( {a{{I( {1,1} )}\mspace{14mu} \ldots \mspace{14mu} {I( {i,j} )}}} )}{L^{{PC}_{g}}( {a{{I( {1,1} )}\mspace{14mu} \ldots \mspace{14mu} {I( {i,j} )}}} )} )}}}$

where L^(PC) ^(g) (a|I(1,1) . . . I(i,j)) is the likelihood of observingthe next pixel with value a conditioned on the previously observed pixelvalues if those pixels belong to grade g tissue and L^(PC) ⁰ (a|I(1,1) .. . I(i,j)) is the true likelihood of observing the next pixel withvalue a.

The smaller the error, the more likely a pixel belongs to a PC class. Itis reasonable, therefore, to weight by the inverse error:

$w_{g,{ij}} = \frac{1}{ɛ_{ij}^{{PC}_{g}}}$

Once the weights are computed for each PC level on each pixel, then thesoft map can be calculated using.

Method for Determining the Most Likely Gleason Distribution

An new method provides a measure of the uncertainty about whether aparticular observed query image pixel is of one pattern grade oranother. This uncertainty is here called a transition probability andrepresents that probability that a pixel is part of one specifiedpattern instead of another specified pattern. This measurement appliedto every pixel within a query image then provides information about theglobal structure of the image—specifically, the uncertainty in thedistribution (or relative population) of pixel grades within a map.Analysis of this uncertainty combined with the observed soft mapdistribution allows for estimation of the most likely Gleasondistribution within the image, which is the amount of each grade patternwithin the image.

Specifically, the percentages of dominant patterns (as determined fromeither hard or soft maps) along with their computed transitionprobabilities generate a prediction of the most probable patterndistribution within the sample.

Secondary likelihoods are given by the weight of one Gleason patternwhen another is dominant. The total transition probability from thedominant pattern to the secondary pattern is the sum of the weights ofthe secondary pattern where the other is dominant divided by the totalsum of all pattern weights where the other is dominant Using the 3, 4,and 5 patterns only, the transition probability T_(d→s) from a dominantpattern d into secondary pattern s is then

${T_{darrow s} = \frac{\sum\limits_{i,j}^{\;}\; {c_{d,{ij}}w_{s,{ij}}}}{\sum\limits_{g,i,j}^{\;}\; {c_{d,{ij}}w_{g,{ij}}}}},{where}$$c_{dij} = \{ {{{\begin{matrix}1 & {w_{d,{ij}} > w_{{g \neq d},{ij}}} \\0 & {{otherwise},}\end{matrix}{for}\mspace{14mu} {each}d} = 3},4,{5;{s = 3}},4,{5;{{{and}g} = 3}},4,5} $

The transition probabilities from each pattern into another can becollected into a single transition matrix:

$T = \begin{bmatrix}T_{3arrow 3} & T_{4arrow 3} & T_{5arrow 3} \\T_{3arrow 4} & T_{4arrow 4} & T_{5arrow 4} \\T_{3arrow 5} & T_{4arrow 4} & T_{5arrow 5}\end{bmatrix}$

Define a density of states vector M(k) which contains the observablemembership percentages of each pattern within a biopsy sample atiteration k. Each iteration describes a more likely distribution thanthe last based on the transition probabilities. If there are n₃ hardgrade 3 pixels, n₄ hard grade 4 pixels, and n₅ hard grade 5 pixels in aquery image, then the image's initial density of states is

${M( {k = 0} )} = {\begin{bmatrix}{\% \mspace{14mu} {of}\mspace{14mu} {Grade}\mspace{14mu} 3\mspace{14mu} {Patterns}} \\{\% \mspace{14mu} {of}\mspace{14mu} {Grade}\mspace{14mu} 4\mspace{14mu} {Patterns}} \\{\% \mspace{14mu} {of}\mspace{14mu} {Grade}\mspace{14mu} 5\mspace{14mu} {Patterns}}\end{bmatrix} = \begin{bmatrix}\frac{n_{3}}{n_{3} + n_{4} + n_{5}} \\\frac{n_{4}}{n_{3} + n_{4} + n_{5}} \\\frac{n_{5}}{n_{3} + n_{4} + n_{5}}\end{bmatrix}}$

A prediction for the density of states at time k is given by

M(k)=T ^(k) M(0)

(if and only if the transition matrix is stationary with respect to k).If all transition probabilities are greater than zero, then T is ergodicand models an aperiodic Markov chain which has a unique, steady statesolution for any initial density of states with nonzero entries:

M(∞)= V _(λ=1)

where V _(λ=1) an eigenvector of T for eigenvalue λ=1 (which in thiscase is guaranteed to exist) with entries that sum to 1. Because theeigenvalue of the eigenvector equals 1, successive operations of thetransition matrix T on V _(λ=1) can only yield V _(λ=1), which must bethe most likely distribution of Gleason patterns within a samplecomparable to the query image. Under these assumptions, the most likelydistribution of Gleason patterns within query sample are completelycharacterized by observed secondary transition probabilitiesirrespective of the initial distribution of patterns observed within thesample.

System for Prostate Cancer Gleason Grading/Scoring Based on Morphologyand Architecture of Prostatic Tissue (System 1)

Referring now to another embodiment of the invention, this systemprovides several methods including methods for tissue objectssegmentation, gland reconstruction, morphology and architecturalfeatures extraction, and image classification. The proposed systemproceeds into two phases: (1) the description (Color, tissue structures,gland and architectural features) that extracts the characteristics fromthe biopsy image pattern and (2) the classification that recognize theGleason grades by using some characteristics derived from the firstphase by implementing learning algorithms. The flow diagram of theproposed system is presented in FIG. 7.

Referring now to this embodiment of the invention in more details, theimage-level features include, color information such as colordistribution entropy and 2-D color histograms. To compute 2-D histogramsfeature vector a color model transformation may be perform. Thistransformation converts each RGB image into its equivalent form in asuitable color space; for instance, L*a*b or HSV (Hue, Saturation andIntensity) color models. Next, 2-D histograms are computed using n binsat each dimension. The number of bins may be selected by variousheuristic methods, for example, cross-validation.

In further detail, still referring to the morphology- andarchitectural-based grading system, tissue structure are segmented usingany suitable segmentation method, for example the color region mappingproposed in another embodiment of this invention. Having segmentedtissue structure and gland units the following features may be extractedfrom the objects: First order statistics of morphology of lumen andgland such as size and shape, gland boundaries features, nuclei densityat several scales, nuclei area and circularity. In addition,architectural features from graphs, for example Voronoi, Delaunay andMinimum Spanning Tree (MST) graphs constructed from color primitives(after segmentation procedures) may be extracted. Such architecturalfeatures include but are not restricted to number of nodes, statisticsand entropy of edge length and polygon area, eccentricity, fractaldimension, among others. Once the feature extraction process isfinished, algorithms for feature selection and ranking may beimplemented to choose the features that best represent the tissue forclassification purposes. The selected features are processed such thatthey are scaled according to their importance.

Referring now to the image classification step, several learningalgorithms may be used. The classifiers may be decision trees, k-nearestneighbor, Bayesian, neural networks, and/or SVM. Boosting algorithmssuch us AdaBoost, SVM Boost, etc. can be used to improve theclassification performance. Next table shows the classification accuracyranges of the multiclass recognition task for grading prostatecarcinomas regions of interest.

Gleason Correct classification rate pattern for different classes ofimages Grade 3 93.33%-99.9% Grade 4 86.67%-99.9% Grade 5 100%

Textural-Based System for Prostate Cancer Gleason Grading and Scoring(System 2)

In this embodiment, transform is considered for transform (wavelets,cosine, Fourier, and others) features extraction for classification ofprostate cancer biopsy images. The flow diagram of this system for awavelet transform is outlined in FIG. 8. Images are preprocessed usingsome of the previously discussed preprocessing techniques in theprevious embodiments prior to image transform and features extraction.Features are classified using any suitable learning algorithm forexample Support Virtual Machine (SVM). Classification performance showedthat Haar wavelet transform produced an improved result in terms ofenergy features extractions when images are first preprocessed using theconcept of morphological filtering and transform coefficients ratio.

Wavelet Transform: In wavelet transform, the signal is decomposed into afamily of orthogonal ψ_(m,n)(x) and ψ(x) as follows:

ψ_(m,n)(x)=2^(−m/2)ψ(2^(−m) x−n) for integers m and n

For the construction of the mother function ψ(x), the determination of ascaling function Φ(x) is needed to satisfy the two-scale differenceequation:

This yields the wavelet kernel ψ(x) as follows:

${{\varphi (x)} = {\sqrt{2}{\sum\limits_{k}^{\;}\; {{h(k)}{\varphi ( {{2\; x} - k} )}}}}};$${{\psi (x)} = {\sqrt{2}{\sum\limits_{k}^{\;}\; {{g(k)}{\varphi ( {{2\; x} - k} )}}}}};$

where g(k)=(−1)^(k) h(1−k)

For wavelet decomposition to a Jth-level, the relation of the decomposedcoefficients to direct previous coefficients can be given as:

c _(j+1,n)=Σ_(k) c _(j,k) h(k−2n) and d _(j+1,n)=Σ_(k) c _(j,k) g(k−2n)where 0≦j≦J

For 2-D wavelet transform, the basis functions operate along thehorizontal and vertical directions as a pair of low pass and high passfilters, h and g, respectively. The coefficients at the HH, HG, GH, andGG correspond to the Approximated, Horizontal, Vertical, and Detailedsubmatrices in the 2D transform space.

To illustrate the use of a transform in this embodiment, the Haarwavelet transform is shown as an example. The Haar basis function can beexpressed as follows:

${{H(2)} = {{\begin{pmatrix}{+ 1} & {+ 1} \\{+ 1} & {- 1}\end{pmatrix}\lbrack{Haar}\rbrack}_{2^{n}} = {{H( 2^{n} )} = \begin{pmatrix}{{H( 2^{n - 1} )} \otimes ( \begin{matrix}{+ 1} &  {+ 1} )\end{matrix} } \\{\sqrt{2^{n - 1}}{{I( 2^{n - 1} )} \otimes ( \begin{matrix}{+ 1} &  {- 1} )\end{matrix} }}\end{pmatrix}}}},{n = 2},3,\ldots \mspace{14mu},$

Many different features can be extracted from wavelets and orthogonaltransforms such as cosine, sine, Hartley, Fourier, and others. Someexemplary features may include the normalized coefficients of energy andentropy features for each decomposed submatrix at a given level asfollows:

${energy} = \frac{\sum\limits_{i}^{\;}\; {\sum\limits_{j}^{\;}\; x_{ij}^{2}}}{nxn}$${entropy} = {{- ( \frac{1}{nxn} )}{\sum\limits_{i}^{\;}\; {\sum\limits_{j}^{\;}{( \frac{x_{ij}^{2}}{{norm}^{2}} ){\log ( \frac{x_{ij}^{2}}{{norm}^{2}} )}}}}}$

Where x_(ij) represents wavelet coefficients; n is the dimension of thesubmatrix, and

${norm}^{2} = {\sum\limits_{i}^{\;}\; {\sum\limits_{j}^{\;}\; x_{ij}^{2}}}$

Due to the diversity of staining colors at each image, the performanceof these features should be evaluated against each color component ofany color model, for example RGB color. Energy and entropy features arethen extracted for each submatrix and classified using the SupportVector Machine (SVM) classifier as an example. Features in the trainingset contain one class label and several multiple features. The learningobjective is to sort the training data into groups/classes so that thedegree of association is strong among the features of the same class andweak among other members of other classes. The set of features in ourcase represent the wavelet-based feature vectors.

As discussed previously, there is significant classification performancediversity among the different grades. In order to improve theclassification performance, the concept of deploying the Human VisualSystem (HVS) rational filters is considered in both, the spatial domainof color channels, and in the wavelet transform domain.

In previous embodiment, we developed an algorithm for visualizing anddetecting edges in a colored image using logarithmic approach. Some ofthe proposed morphological filters include the following:

${\beta \; 1( \frac{\max ({Fn})}{\max ({Fm})} )};$${\beta \; 2( \frac{\min ({Fn})}{\max ({Fm})} )};$${\beta \; 3( \frac{\min ({Fn})}{\min ({Fm})} )};$$\beta \; 4( \frac{\max ({Fn})}{\min ({Fm})} )$

Where βn is a constant used for grayscale normalization.

For the class of biopsy images, taking as an example the RGB colormodel, the filter considers the logarithmic ratio of the minimum RGBcomponent to the maximum one. This type of filters has the ability tosegregate the dark basal cells from other tissue characteristics. Itoperates on biopsy images as follows:

$I_{new} = {\beta \; {\log_{10}( \frac{\min ({RGB})}{\max ({RGB})} )}^{\alpha}}$

Where I_(new) is the output processed/filtered image.

In another embodiment, the concept of rational filters in the transformdomain is used. The wavelet transform provides a good localizationproperty in space-frequency domain. Since the most significantcoefficients information exist mainly in the approximated and detailedcoefficient sub matrices, the ratio of these two matrices is taken, thenthe energy of the coefficients ratios is considered in classification.

A transform approach for image quality measurement using the propertiesof HVS has been proposed. For any wavelet transform, the rational filtercoefficients for the low-low “Approximated” and high-high “Detailed”frequency filter banks in a 1-D signal can be given as follows:

$C_{{ratio}_{1\; D}} = \frac{\sum\limits_{i}^{\;}\; {Ci}_{A}}{\sum\limits_{i}^{\;}\; {Ci}_{D}}$

Similarly, the coefficients ratio, C_(ratio), for the 2-D case iscalculated as follows:

$C_{{ratio}_{2\; D}} = \frac{\sum\limits_{i}^{\;}\; {\sum\limits_{j}^{\;}\; {Cij}_{AA}}}{\sum\limits_{i}^{\;}\; {\sum\limits_{j}^{\;}\; {Cij}_{DD}}}$

Another proposed method is to use the phase information of the Lo and Hicoefficients ratio, which is expressed for 1-D signal as:

$\theta_{i} = {\tan^{- 1}( {\frac{\sum\limits_{i}^{\;}\; {C_{i}{Lo}}}{\sum\limits_{i}^{\;}\; {C_{i}{Hi}}}} )}$

Features extracted from the proposed method are fused with the waveletenergy feature. The following table tabulates the classificationperformance for different classes of prostatic carcinomas by employingthe explained approach.

Gleason Correct classification rate pattern for different classes ofimages Grade 3 77.0%-99.9% Grade 4 85.3%-99.9% Grade 5  96.0-100%

The classification performance using the proposed method succeeded inoutperforming the previous method where energy feature is extracted fromsome color channels of the biopsy images directly.

Another extracted feature is related to the fractal dimension (FD)measure. One embodiment of the invention provides two developedalgorithms in fractal analysis. The developed algorithms are based onedge detection information as a preprocessing stage as well as quad-treedecomposition algorithms. The classification performance resulting fromadopting such pre-processing algorithm showed an improved performance insome Gleason grades over the classical implementation of the boxcounting method.

Methods for Computation of Fractal Dimension by Improved Box-CountingAlgorithm

The architecture and geometrical features of biopsy images of prostatepresents a rich domain for fractal analysis through the concept offractal dimension (FD). The fractal dimension provides a geometricaldescription of an image by measuring the degree of complexity and howmuch space a given set occupies near to each of its points. The growthof cancer through the development of different architectures associatedwith different grades shows distinct features that can be used forclassification using fractal dimension. Malignant cancerous cellsdemonstrate darker intensity gray level values compared to thesurrounding tissue parts. In this part, fractal features throughDifferential Box Counting (DBC) will be used.

Given a set S in Euclidean n-space, the self similarity of S is obtainedif Nr distinct copies of S exist by r ratio. The following equationprovides the fractal dimension measure “D” as follows:

$D = \frac{\log ( N_{r} )}{\log ( {1/r} )}$

To find the box-counting dimension, the image is binarized by using anyproper approach in order to have two levels of intensities, either 0 or1, rather than 256 gray levels intensities. To find a suitable binaryimage for fractal dimension computation several algorithms such as theproposed logarithm edge detection, gray level thresholding, adaptivethresholding, and bit-plane decomposition may be used. Then, the binaryimage is divided into grids of size r and the number of grid boxes Nrexhibiting self similarities is counted. The next step is to reduce thegrid size, r2 and find the corresponding number of boxes, Nr2. Gridsizes are reduced in a pattern following 2n. For instance, a 512×512image yields the following grid sizes, 512, 256, 128, 64, 32, 16, 8, 4,and 2 grid sizes. The resulting fractal dimensions for each grid sizeare fused together to obtain a unique fractal dimension.

Also, in this embodiment we proposed an edge detection-based fractaldimension measure (EDFDE). It is utilizing the concept of edge detectionusing parameterized logarithmic image processing, which was developed inprevious embodiments. Biopsy images go through PLIP-based edge detectionprior to box counting.

In summary, feature vectors were extracted from wavelet and fractalanalysis domain. Each extracted feature vector represents a differentarchitectural property used for classification purposes feature vectorsfrom wavelet transform provided a high performance for classifying andgrading biopsy images. Before classification step, feature selection isperformed. The objective is to select features sets that best describesthe statistical property of the biopsy images classification features.There are various algorithms for features selection and validation offeatures prior to classification in literature. For example, theminimum-Redundancy Maximum-Relevance (mRMR) features selection methodmay be used. Using the mRMR method, the features space (N-features) isreduced to M-features (N>M).

Referring in more details to the image classification step, severallearning algorithms may be used. The classifiers may be decision trees,k-nearest neighbor, Bayesian, neural networks, and/or SVM. Boostingalgorithms such us AdaBoost, SVM Boost, etc. can be used to improve theclassification performance. The table below shows an example ofclassification of images from different Gleason grades and highlightsthe final SVM classification performance using “leave-one-out”cross-validation.

Gleason Correct classification rate pattern for different classes ofimages Grade 3 93.0%-99.9% Grade 4 89.3%-99.9% Grade 5  96.0-100%

Textural-Based System for Prostate Cancer Gleason Grading and Scoring(System 3)

In this embodiment, the invention provides a system with several methodsfor automatic classification of prostate cancer carcinomas by usingwavelet-based and fractals-based textural features. Image-level featuresare extracted from wavelet decomposition and color fractal analysis ofimages, as well as object-level fractal dimension features as it isshown in FIG. 9.

Referring now to the embodiment in more detail, the following featuresare extracted from the cancerous images:

Transform-Based Features:

In order to compute the feature vectors, real/complex wavelets ormulti-wavelets or other orthogonal transforms (e.g. Fourier, cosine,sine, Hartley, Hadamard) may be used in general. Several wavelettransforms are widely used for quantitative image processing since theyprovide an effective tool for multi-resolution analysis. A particularapproach may include parameterized Slant-Haar transform [19], which isdescribed below:

Slant-Haar wavelet is the wavelet transform of choice for this example.The slant-Haar transform of order 2^(n) for n=3, 4, 5 . . . , is definedby:

X=S ₂ _(n) x

Where x is an input data vector of size n.

S_(2^(n)) = S_(2^(n))(β₄, …  , β_(2^(n)))$S_{2^{n}} = {Q_{2^{n}}\lbrack \frac{A_{2} \otimes S_{2,2^{n - 1}}}{I_{2} \otimes S_{{2^{n} - 2},2^{n - 1}}} \rbrack}$

For −2^(2n-2)≦β₂ _(n) ≦2^(2n-2) and n≧3

Where S_(2,2) _(n-1) is a matrix of dimension 2×2^(n-1) comprised of thefirst two rows of S₂ _(n-1) and S₂ _(n-1) _(-2,2) _(n-1) is a matrix ofdimension (2^(n-1)−2)×(2^(n-1)) comprised of the third to the (2^(n-1))th rows S₂ _(n-1) . The symbol

denotes the operator of the Kronecker product. The matrix A₂ is givenby:

$A_{2} = {\frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1 \\1 & {- 1}\end{bmatrix}}$

And the matrix S₂ _(n-1) is the recursion kernel defined by:

$Q_{2^{n}} = \begin{bmatrix}1 & \mspace{11mu} & \; & \; & \; & \; \\\; & b_{2^{n}} & a_{2^{n}} & \; & \; & \; \\\; & a_{2^{n}} & {- b_{2^{n}}} & \; & \; & \; \\\; & \; & \; & 1 & \; & \; \\\; & \; & \; & \; & 0 & \; \\\; & \; & \; & \; & \; & 1\end{bmatrix}$${a_{2^{n}} = \sqrt{\frac{3( 2^{{2\; n} - 2} )}{{4( 2^{{2\; n} - 2} )} - \beta_{2^{n}}}}},{b_{2^{n}} = \sqrt{\frac{( 2^{{2\; n} - 2} ) - \beta_{2^{n}}}{{4( 2^{{2\; n} - 2} )} - \beta_{2^{n}}}}}$

For −2^(2n-2)≦β₂ _(n) ≦2^(2n-2) and n≧3

Another approach may include parameterized PLIP Slant-Haar transform. InPLIP Slant-Haar wavelet, the arithmetic operations in waveletdefinitions are replaced by PLIP operations [14]:

$ {g_{1} + g_{2}}arrow{g_{1} \oplus g_{2}}  = {g_{1} + g_{2} - \frac{g_{2}g_{2}}{\gamma (M)}}$$ {g_{1} - g_{2}}arrow{g_{1}\Theta \; g_{2}}  = {{k(M)}\frac{g_{1} - g_{2}}{{k(M)} - g_{2} + ɛ}}$g₁g₂ → g₁ * g₂ = ϕ⁻¹(ϕ(g₁)ϕ(g₂))

Where k(M), γ(M) are arbitrary linear functions of the type γ(M)=AM+Bfor constants A and B. The parameter ε is a very small constant includedto avoid division by zero in the case k(M)=g₂.

Referring now in more detail to the features which can be extracted fromwavelet coefficients, the following features are considered:

-   -   Statistics of wavelet detail coefficients: Using a wavelet        transformation of the three color components of the image in one        or more color models at a specific scale, the first order        statistics of scaled coefficients of each sub-band are computed        and arranged as a feature vector. Those statistics may include        average, median, standard deviation, high order moments, among        others.    -   Phase information of approximation and detail coefficients:        After transforming the image using any suitable wavelet basis,        the phase information generate by the approximation and detail        coefficients may be computed as follows:

$\theta_{i} = {\tan^{- 1}( {\frac{\sum\limits_{i}^{\;}\; {C_{i}{Lo}}}{\sum\limits_{i}^{\;}\; {C_{i}{Hi}}}} )}$

-   -   Joint Probability distribution of wavelet detail coefficients:        The color input image is first decomposed into its RGB        components. Before applying 2-D wavelet filters, color model        transformation can be performed. Next, wavelet decomposition is        performed for each channel separately and joint probability of        scaled detail coefficients is computed for each possible pair of        channels. The joint distribution consists of a set of bins which        counts the number of wavelet details coefficients at a given        scale falling in a given area of the plane:

${p_{{Chi},{Chj}}( {x,y} )} = {k{\sum\limits_{x}^{\;}\; {\sum\limits_{y}^{\;}\mspace{14mu} {{AND}\mspace{14mu} ( {{{Ch}_{i} = x},{{Ch}_{j} = y}} )}}}}$∀i ≠ j; i, j = 1, 2, 3

In equation 3, i,j=1, 2, 3 indicate the three channels of a color inputimage. x, y represent the numerical value of the coefficients whichrange, in our experiments, from 0 to 255. The logical AND function willreturn 1, when both arguments are true, and 0 otherwise in order tocount the number of coefficient falling in a given bin.

-   -   Wavelet Generalized Energy Distribution: In order to capture how        the wavelet energy is distributed within the images, each input        image (gray scale image or R-G-B components) is divided into        non-overlapping s×s blocks. A wavelet transformation, for        example Haar transform, is then applied to each block        independently. The wavelet transform converts the data array        into a series of wavelet coefficients, each of which represents        the amplitude of the wavelet function at a particular location        within the array and for a given wavelet scale. The generalized        energy of approximation sub-band can be computed as:

${E_{i}^{A} = {\sum\limits_{j = 1}^{N}\; {A_{ij}}^{p}}},{i = 1},\ldots \mspace{14mu},K,{P = 1},2.3,\ldots$

Here, K is the wavelet decomposition level. N is the number of thecoefficients in the approximation sub-band at each decomposition level.First order statistics, i.e. mean, standard deviation, median, minimumto maximum ratio, etc, of the energy distribution are computed toproduce feature values.

According to an embodiment of this invention, the sub-band coefficientsare processed by a window function, wherein the center coefficientx_(ij) is replaced by the result of the function ƒ_(w)(x_(i,j→m×n)) in asmall neighbor w of size m×n, m,n=1, 2, 3, . . . .

For example, the function for a 3×3 neighbor could be:

ƒ_(w)(x _(i,j→3×3))=Σ[x _(i,j)(x _(i−1,j−1) ^(α) ¹ x _(i−1,j) ^(α) ² x_(i−1,j+1) ^(α) ³ x _(i,j+x) _(i,j) ^(α) ⁴ x _(i,j) ^(α) ⁵ x _(i,j−1)^(α) ⁶ x _(i+1,j−1) ^(α) ⁷ x _(i+1,j) ^(α) ⁸ x _(i+1,j+1) ^(α) ⁹ )]²

Where α_(i) are user-defined constants, for example 0, 1, 2, . . . .

Or,

f_(w)(x_(i, j → 1 × 1)) = ∑x_(i, j)², f_(w)(x_(i, j → 1 × 1)) = ∑Log²x_(i, j), f_(w)(x_(i, j → 1 × 1)) = ∑x_(i, j)Logx_(i, j)f_(w)(x_(i, j → 1 × 2)) = ∑(x_(i, j)x_(i, j − 1))²$\theta_{i} = {\tan^{- 1}( {\frac{\sum\limits_{i}^{\;}\; {C_{i}{\log ({Lo})}}}{\sum\limits_{i}^{\;}\; {C_{i}{\log ({Hi})}}}} )}$$\theta_{i} = {\tan^{- 1}( {\frac{\sum\limits_{i}^{\;}\; {C_{i}{\log ( {{Lo}( {x_{i,j}x_{i,{j - 1}}} )} )}}}{\sum\limits_{i}^{\;}\; {C_{i}{\log ( {{Hi}( {x_{i,j}x_{i,{j - 1}}} )} )}}}} )}$

The features above can be computed using the transformed waveletcoefficients. In the case complex wavelet are used, x_(ij) refers to themagnitude of the coefficients.

Method for Computation of Color Fractal Dimension

In accordance with an embodiment of this invention, a new algorithm forcomputing fractal dimension of color images is proposed. The modifiedalgorithm is based on a color space transformation of the histopathologyimages. To this end, a new modification in the (R,G,B) color space isused. The proposed algorithm extracts a single fractal dimension of acolor image or a fractal dimension vector corresponding to the colorfractal dimension of non-overlapping blocks of size m×n. In thisalgorithm the color image is considered a set of points in a5-dimensional space S. Each point sεS is represented by its two spatialdimensions and its three color components in the new color models_(i)=(x,y, ƒ₁(R,G,B,α₁), ƒ₂(R,G,B,α₂), ƒ₃(R,G,B,α₃)) where ƒ_(i) is afunction of 3 colors and some set of constants α_(i). The functionshould be adjusted to enhance the separation between cancer grades.Various examples of the color transformation are shown in FIG. 10. Colorfractal dimension D of an image can be approximate by:

${N(L)} = {{g(L)}{\sum\limits_{m = 1}^{N}\; {\frac{1}{m}{P( {m,L} )}{\infty L}^{- D}}}}$

Where, N(L) is the total number of boxes (size L) needed to cover theimage, and it is multiplied by a function g(L), such as g(L)=1, L, L²,log L, L log L. P(m,L) is the probability matrix having m pointsincluded into a hyper-cube (or box) of size L centered at an arbitrarypoint of S. Below, the probability matrix calculation procedure isoutlined:

1. Count the number of points sεS,

s _(i)=(x,y,ƒ ₁(R,G,B,α ₁),ƒ₂(R,G,B,α ₂),ƒ₃(R,G,B,α ₃)) with

d (F,Fc)≦L in which d (F,Fc) is a distance between the center pixel Fcand the rest of pixels F in the box, and ƒ_(i), ƒ_(ci), is the value ofthe i-th dimension of point sεS in the 5-D space s=(x,y,r′,g′,b′). Forexample, the distance d (F,Fc) could be:

The Minkowski distance

d(F,Fc)=|F−F _(c)|=max(|ƒ_(i)−ƒ_(ci))∀i=1,2 . . . ,5

Or, p-norm distance

${d( {F,{Fc}} )} = {{{F - F_{c}}} = ( {\sum\limits_{i = 1}^{5}\; {{f_{i} - f_{ci}}}^{p}} )^{1/p}}$

2. For an arbitrary size L normalize the matrix P(m,L) by using thefollowing relationship

${\sum\limits_{m = 1}^{N}\; {P( {m,L} )}} = 1$

in which N is the number of pixels included in a box of size L.

Finally, the fractal dimension D is found by the linear fitting of thenegative of the logarithm of L against the logarithm of itscorresponding. The best function g(L), can be found using

$G = {\max\limits_{g}\mspace{14mu} \{ {\min\limits_{i,j}( d_{i,j} )} \}}$

Where d_(i,j) is the distance between pairwise averages Aν_(i) with i=3,4, 5, of the fractal dimension of the set of class of images of Gleasongrades 3, 4 and 5 calculated for each function g(L).

Simplifications of the invented algorithm result in the color fractaldimension introduced by Ivanovici 20].

Referring now to the image classification step, several learningalgorithms may be used. The classifiers may be decision trees, k-nearestneighbor, Bayesian, neural networks, and/or SVM. Boosting algorithmssuch us AdaBoost, SVM Boost, etc. can be used to improve theclassification performance. Next table shows the classification accuracyranges of the multiclass recognition task for grading prostatecarcinomas regions of interest.

Gleason Correct classification rate pattern for different classes ofimages Grade 3  99.9%-100% Grade 4 86.67%-99.9% Grade 5 100%

3-D Gleason Grading System (System 4)

In this embodiment, 3-D image is obtained using a single 2-D tissueimage to assign Gleason grade and score to the third dimension of thereconstructed image. The building blocks of this part are shown in FIG.11 and will consider the following steps:

(1) Mapping and pre-processing: feature vectors algorithms related towavelet transform, and fractal analysis will be considered. The initialstep considers mapping 2-D images into slices of layers usingpre-developed mapping algorithms in order to have a 3-D voxelrepresentation of a single 2-D image. Preprocessing algorithms willconsider: Haar-HVS based algorithms, logarithmic RGB modelrepresentation, and edge detection-based algorithms.

(2) Feature vectors normalization: feature vectors generated from abovealgorithms will be normalized.

(3) Feature Selection and Classification: feature selection algorithmsmay be deployed to select dominant features generated as describedabove. Several classifiers can be used for recognition. For example SVMclassifier may be used to classify such features and generate a Gleasonscore.

(4) Integration and validation: the final output of this subsystem maybe verified and integrated with the 2-D Gleason subsystem.

Referring now in more details to this embodiment of the invention, themapping of a single 2-D image into a 3-D space may be done by usingseveral mapping algorithms. The used algorithm are grouped in depth cuessuch as defocus, linear perspective, shading, patterned textures,symmetric patterns, occlusion, statistical patterns, among others.

System for Whole-Slide Processing and Prostate Cancer Diagnosis:Detection, Grading, Mapping, Scoring (System 5)

A computer-aided diagnosis system addresses the detection, Gleasongrading, scoring and mapping of prostate carcinoma in whole digitizedslides. The system follows a multi-resolution approach, in whichspecific features are extracted at a specific image magnification. Forexample: (1) at low resolution color and textural features areextracted, (2) at medium resolution, architectural tissue analysis isperformed via graph features, and (3) At high resolution, themorphological and wavelet features are used for tissue characterization.The system proceeds as follows: The biopsy specimen is segmented fromthe background to save calculation time, a square or rectangular grid ofsize m×n is superimposed for patch-based cancer diagnosis (the size ofthe grid could change depending on the resolution), spatial andtransform domain (multi-resolution analysis) features are extracted andblock classification is performed by using any of the systems of theprevious embodiments.

A classification system according to an embodiment of the invention usesmultiple features for classification of Whole-Slide and several computeralgorithms which produce the following final outputs: tumor extent inthe needle biopsy, number and identification of positive needle cores,percentage of positive cores for cancer, length of the tumor in positivecores and length of the core, percentage of needle core tissue affectedby carcinoma and area of tissue positive for cancer, Gleason grade ofeach selected cancerous region of interest, percentage of each Gleasonpattern in biopsy specimens, Gleason score of each positive core,including tertiary pattern information, the most frequent Gleason scorein slides, numbers of identical Gleason scores in the positive cores,and localized cancer map including the grade of each cancerous patch inall needle biopsy slides as shown in FIG. 12.

In accordance to this embodiment, the classification algorithm usesmorphology-, textural-based classification subsystems or a combinationof those systems provided in previous embodiments of the presentinvention. As a result, a grade of the Gleason grade scale is assignedto each cancerous patch in the whole slide and the Gleason score iscomputed automatically. In addition, the system is able to detecttertiary patterns present in the biopsy image. The percentage of extentof each grade found in the image is also computed.

The system of studying biopsy samples may be an automated procedure toprocess large number of biopsy slides by integrated a production linemodel. The automated system output automatically classifies 2-D and 3-DGleason grades by using serial sections/slices of prostate core biopsyimages. The automated system output automatically recognize de Gleasongrade, compute the Gleason score and predicts the likelihood of prostatecancer reoccurrence after treatment, as well as more accurate riskassessment, and the decision on the necessity of further biopsy samples.The automated system output allows improved and efficient treatmentplanning.

In an embodiment, the method and system allow a surgeon to beimmediately aware of whether the surrounding tissue of a specimencompletely surrounds the tumor, prior to the time that the patient isclosed up.

3-D Core Grading/Scoring System (System 6)

Digital 3-D reconstruction of tissue has many uses in medical researchand it can be also used for cancer diagnosis because more features canbe analyzed from high-resolution 3-D images. For example, additionalinformation regarding the volume of tumor areas, the distribution ofcells, the branching of blood vessels in a tumor, etc can be extractedto make better and more accurate cancer diagnosis and prognosis.

The building blocks of this embodiment will consider the following:

(1) Acquiring images: serial sets of images (M sets from N patients with12 slices each) need to be pre-labeled by a multiple experiencedpathologists with Gleason score ranging from 3+3 to 5+5 for learningalgorithm training.

(2) Building 3-D model: building a 3-D image from a set of slices forpreprocessing and featuring vectors generation in the next step.

(3) Preprocessing and feature vectors generation: preprocessingalgorithms presented earlier in this proposal are used for theinvestigation of the core pattern and a possible feature vectorsgeneration. The extracted features may include: shape and structure ofcells, graph features, ‘micro-architecture’ of blood vessels, andtumors' volume, shape, regularity, among others.

(4) Classification and validation: in this part, a Gleason grade andscore are given to the 3rd dimension of the core based on SVM classifieror other suitable classifiers with different validation methods. Thelabeled grades are investigated and verified with a pathologist toinsure validation of the grade and its correlation with the 2-D gradingsubsystem.

In accordance to an embodiment of this invention, the image acquisitionprocess may be performed by using a photon-light, a stereo-lightmicroscope or an electron microscope. The selection of the microscopeshould be done according to the required image resolution.

To construct 3-D image from 2-D slices, different software tools can beused. For instance, ImageJ, Fiji, or developed algorithms in Maltab,among others. The general flow diagram for 3-D reconstruction is shownin FIG. 13.

According to another embodiment of the invention, tissue imagepreprocessing and structure segmentation may be performed using 2-Dslices. FIG. 14 shows an example of 3-D reconstruction of segmentedlumens and epithelial nuclei from a prostate cancer biopsy image.Feature extraction may be performed from both 2-D and 3-D images.

According to another embodiment of the present invention, Gleasongrading and scoring is performed in the 3-D tissue image fromlongitudinal and transverse views (see FIG. 15). The total volume of thetumor is also calculated. In this part, the Gleason score of the tumorwill be computed as the combination of the most predominant grades inboth directions. For example, the most predominant grade could be theaverage or mode of the grades which were found in the longitudinal andtransverse tissue evaluation. If the distribution is bimodal ormultimodal, the first component of the 3-D Gleason grade will be themaximum grade among the modal values and the second component will bethe next grade value. Ternary patterns may be included in the finalGleason score, especially if those patterns correspond to a high Gleasongrade.

System for Whole-Slide Processing and Prostate Cancer Diagnosis:Detection, Grading, Mapping, Scoring (System 7)

A computer-aided diagnosis system addresses the detection, Gleasonscoring and mapping of prostate carcinoma in whole digitized slides. Thesystem follows a prediction/compression (PC) approach, in whichpredictors are trained to optimally compress different Gleason patterns.For example: (1) segment the query image (color, tissue structures,gland and architectural features) into pattern-defining characteristics,(2) reduce query to quarter scale using a wavelet transformapproximation, (3) unravel 2-D query into 8 distinct space-fillingcurves, (4) train a universal predictor/compressor (UDC) on each Hilbertcurve, (5) train 3 UDCs on collections of already-graded 3, 4, and 5patterned, archival images, (6) compress the query image Hilbert curveswith each UDC, (7) statistically analyze compression results, and (8)grade query image and pixels with Hard and Soft scores computed from thestatistical results.

Biopsy image preparation is necessary for formatting image data into astructure usable by the novel universal compressor-detectors (UDCs).Images need to be simplified and reduced to support fast data analysisand also need to be one-dimensionalized for construction andimplementation of a UDC. Specifically, image data must support twophases of UDC action:

1. Training and

2. Detection.

Training builds UDC structures on pathologist-graded,one-dimensionalized, Gleason pattern images. Detection involvesprediction and compression operations by a UDC on a similarly preparedquery image.

Thus, image preparation for both training images and query imagesconsists of three significant steps: (1) structural segmentation forimage simplification, (2) image downscaling for data reduction, and (3)1D conversion for UDC support. While different approaches may be takento perform these steps, the methods chosen here are designed tofacilitate computational speed while simultaneously maintaining therobustness of the system.

First, 512×512 pixel, H&E-dyed biopsy images are segmented into sixsignificant structural elements including: lumen, epithelial cytoplasm,stroma, blue mucin, and epithelial nuclei. The method simultaneouslyextracts the important features found in the various Gleason patternsand reduces the dynamic range of pixel values (originally 0-255 in eachof the red, green, and blue channels) to 0-5, effectively reducing theimage alphabet size from the approximately 1.68×10⁷ color combinationsof an RGB image to 6, which is a very practically sized alphabet for aUDC to handle.

Images are also reduced to quarter scale (64×64 pixels) in order toreduce the number of pixels from 262,144 per image to 4,096 per image.This action preserves significant regional statistics within the sampleimages while alleviating the amount of data on which a UDC must train orcompress. Scale reduction is achieved using a 2-level integer Haarwavelet transform to generate a 64×64 approximation of the original512×512 segmented images and to preserve the dynamic range of thealphabet. Other wavelets could be used for the same task; however, theHaar is chosen because of its ability to be computed with extreme speed.

Lastly, the scaled, segmented image is unwrapped into one dimension tosupport use with a UDC. For example, this can be achieved by usingseveral Hilbert curve raster scans. A Hilbert curve is afractally-generated, space-filling curve on integer powers of 2 sidedhyper-squares. For example, for every square digital image with sides of2′ (where p is an integer) pixels a Hilbert curve 25] can be constructedwhich traverses all pixels continuously through adjacent pixels. Infact, 8 unique Hilbert curves can be constructed for square 2-D images.

For UDC training (and compression), a 64×64 Haar approximation of animage is unwrapped into its 8 unique Hilbert curves on which UDCs can betrained. For query image compression, a UDC compresses each of the 8Hilbert sequences separately, and the average compression rates for eachpixel are recorded. Thus, if UDC_(g) predicts the probability p_(ijk)^(g) for pixel (1,1) from the context of kth Hilbert curve (k=1, 2 . . .8), then the observed compressed size in bits of pixel (i, j) usingUDC_(g) is

$Y_{ij}^{{UDC}_{g}} = \frac{\sum\limits_{k = 1}^{8}\; {- {\log_{2}( p_{ijk}^{g} )}}}{8}$

The distributions of average pixel compression rates for a query imageusing each UDC are used to determine the Gleason grade of the queryimage. In the following sections, reference to a query image or itspixels specifically refers to the segmented, down-scaled versions andnot the original query image.

The novel grading system supplies both hard grades (grades where thedecision must be either grades 3, 4, or 5 and no in between) and softgrades (where the grade is a fuzzy membership average of multipleclasses, generating a continuous scale Gleason grade) to a query imageas a whole based on UDC compression rates. Query image pixels are alsoassigned individual hard and soft grades depending on the “distance” ofa trained UDC's compression model to an optimal compression model afterobservation of that pixel. (Note that only images of grades 3, 4, and 5are considered as these represent the significant patterns to cancergrading and were available for our analysis.)

The system roughly consists of roughly 8 discrete parts, illustrated inFIG. 9:

1. UDC Training: 3 UDCs are trained on respective classes of images. 1UDC is trained on a query image under inspection.

a. UDC₃ is trained on pathologist-graded pure pattern 3 images.

b. UDC₄ is trained on pathologist-graded pure pattern 4 images.

c. UDC₅ is trained on pathologist-graded pure pattern 5 images.

d. UDC₀ is trained on a single query image.

2. Query Image Compression: Each UDC compresses an ungraded query image,I, and each pixel (i,j)'s compressed size in bits is stored as anobservation Y_(ij) ^(UDC) ^(g) . Concurrently, pixel errors ε_(ij)^(UDC) ^(g) are also computed using the Kullback-Leibler divergence of apattern UDC with respect to UDC₀, which is the optimal compressor forthe query image by construction.

a. UDC₃ compresses I with compressed pixel sizes Y_(ij) ^(UDC) ³ .

b. UDC₄ compresses I with compressed pixel sizes Y_(ij) ^(UDC) ⁴ .

c. UDC₅ compresses I with compressed pixel sizes Y_(ij) ^(UDC) ⁵ .

d. UDC₀ compresses I with compressed pixel sizes Y_(ij) ^(UDC) ⁰ .

e. KL divergence measures the model errors ε_(ij) ^(UDC) ³ of UDC₃ toUDC₀ produced by observation of pixel I_(ij).

f. KL divergence measures the model errors ε_(ij) ^(UDC) ⁴ of UDC₄ toUDC₀ produced by observation of pixel I_(ij).

g. KL divergence measures the model errors ε_(ij) ^(UDC) ⁵ of UDC₅ toUDC₀ produced by observation of pixel I_(ij).

3. Analysis of Variance Testing: ANOVA is used to determine whether ornot the mean compression rates provided by UDC₃, UDC₄, and UDC₅ arestatistically dissimilar.

4. Hard Grading: A hard decision H for the query image Gleason grade ismade corresponding to the best-compressing UDC's class.

5. Hard Mapping: The query image is mapped, pixel by pixel, with a hardgrades H_(ij) corresponding to the minimally divergent UDC model at eachpixel.

6. Soft Grading:

a. If ANOVA testing indicates the UDC compression rates are notstatistically dissimilar, the query image is given a soft gradeS=UNLCASSIFIABLE, meaning the image is not a distinct member of any UDCclass.

b. If ANOVA testing indicates that at least 2 UDCs are statisticallydissimilar

i. Tukey's multiple comparison test is used to test the independence ofeach UDC's mean compared to the others' means.

ii. The query image is given a soft grade S using weighted means basedon the intersections of the t-distributions modeling the mean UDCcompression rates.

7. Soft Mapping: The query image is mapped, pixel by pixel, with a softgrades S_(ij) computed using weighted means based on the inverse KLdivergence of each pixel.

8. Prediction of Most Likely Pattern Composition:

a. The probability of grade transitions from each grade to the othersare computed from the errors associated with each UDC model.

b. The final pattern composition (in percent of each pattern) iscomputed through steady state analysis of the transition model.

Prognosis System for Integration of Biopsy Classifier Systems withPatient and PSA Information (System 8)

In one embodiment of the present invention, a computerized system isdeveloped by introducing a unique recognition and prognostic system thatprovides the following: (1) Recognition and classification of Gleasongrading in 2-D biopsy images, 3-D biopsy images, and 3-D core patterns;(2) faster and accurate classification with the introduction of imagepreprocessing algorithms; and (3) integration of patient personal dataand prostate-specific antigen (PSA) information for more accuraterecognition and prognostics analysis.

Described herein is a system for the recognition and classification ofGleason grades of pathological biopsy images using three sub-systems forhistology grading that are integrated at a later stage with patientinformation such as PSA information, age, race, and family history. Thegeneral hierarchy of such integration is illustrated in FIG. 17.

The first histology subsystem considers the classification of 2-D biopsyimages; one or more of the previous systems may be used. The secondsubsystem, considers the classification of 3-D biopsy images. In thispart, a 2-D image is mapped using certain mapping algorithms to 3-Dspace to have a 3rd dimension. The conventional way of visualizingregular 2-D biopsy images is replaced by a 3-D view to observehistological details in 3-D space. The objective of this subsystem is toverify the classification of the 2-D subsystem in terms of analyzing theunforeseen information that may exist in the 3-D representation. Theclassification performance of the 2-D and 3-D subsystem is going to befused for accurate recognition performance. The third subsystemintroduces the concept of 3-D core grading of prostate biopsy specimen.In this subsystem serial images of a stack of multiple 2-D tissue biopsyslices are captured from a given suspicious lesion in the prostategland. The 3rd dimension of this stack represents the depth of the tumoror “core”. It can be visualized as a 2-D internal image of the prostategland when a vertical cut is made on the stack of the 3-D slices. Theobjective of this subsystem is to investigate tertiary patterns, andassigns a malignancy Gleason grade for the core (depth of the specimen).The final Gleason score resulted from these three subsystems can belooked at as a Gleason grade composed from 3 scores, the first twoscores represent the Gleason score from the first two subsystems, andthe third score is the core score as follows:

Gleason score=(predominant pattern grade+second predominant patterngrade+3-D core grade)

Regarding the PSA and patient information, information from these blockswill be integrated as an additional feature for accurate classificationand prognostics purposes.

Biopsy images of prostate cancer form diversities in terms of structure,color information, edge information, and visualization of fine detailsassociated with each distinct grade. A typical image recognition systemfor classification purpose is composed of the following three mainblocks/phases: (1) Images preprocessing; (2) features extraction(feature vectors); and (3) classification. Considering the first blockthat is concerned with image preprocessing, image details are collectedfrom images at different Gleason grades prior to feature vectorsextraction phase. Preprocessing of images provides significantinformation in terms of visualization, segmentation, edge informationprior to categorizing images in different cancer grades through featurevectors generation and classification.

In one embodiment of the present invention, methods for 2-D and 3-Dgrading systems integration are provided. The steps to achieve aprostate cancer prognosis system as a follows:

(1) Acquire PSA information attached with each graded biopsy image asfollows: PSA level, PSA-Velocity, % free-to-total PSA.

(2) Acquire patient information attached with each labeled biopsy imageas follows: age, race, Body-mass Index (BMI), and family history ofcancer.

(3) Feature Vector Generation and Classifications of the aboveinformation for features integration with biopsy images featuresassociated with 2-D, 3-D and Core grading classification. The systemsprovides in previous embodiments may be used.

(4) Verification of the whole system using subjective and/or objectivemeasures. System output may be compared with biopsies graded by expertpathologists. Mean Opinion Score (MOS) may be used to measure theperformance of the system. Also, DNA Microarrays results, if availablefor each patient, may be analyzed and compared with our system finaloutput as an objective measure.

The introduction of a new novel Gleason grading concept consideringdifferent aspects of biopsy images ranging from 2-D to 3-D contributesto an enhanced output of more than just an accurate classification ofcancer grade. The system of this embodiment integrates PSA and patientinformation along with results of cancer automatic grading from theproposed histology system. Hence, the final “expert system” output isexpected to provide the following functions:

-   -   Provides an automated Gleason score to assets pathologist in        accurate classification of cancer grades of prostate biopsy        images.    -   Aids medical doctors by providing prognostic information about        the stage and level of prostate cancer.    -   The generated accurate Gleason score contributes directly to        more accurate risk assessment readings using the common        available risk assessment calculators. Accuracy Contributes to        more accurate risk assessment by integrating a more accurate        Gleason grade with other risk assessment variables.    -   Biopsy planning by predicting biopsy locations having suspicious        malignant tissues.

Implementation of the method/or system of embodiments of the inventioncan involve one/or combination of pre-developed or new risk assessmentmodels and 2-D and 3-D Gleason grading/scoring systems. For example, wemay use one of following risk assessment model: Prostate CancerPrevention Trial Risk Calculator (PCPTRC) 21], D'Amico model 22], Kaftannomograms 23], and Cancer of the Prostate Risk Assessment (CAPRA) 24]).In general, a risk assessment system uses PSA level (blood test),Gleason grade, and T stage (size of the tumor on rectal exam and/orultrasound) to group men as low, intermediate, or high-risk. In allthese methods, the accurate detection, grading and scoring of biopsyimages is a crucial step towards a better risk assessment.

System and Methods for Estimation, Prediction, and Cancer GradingOutside of the Histopathology Biopsy Images

In another aspect of this embodiment, a method is provided forpredicting the Gleason grade of blocks in histopathology biopsy images.The basic idea is to predict the state of missing blocks or to produceextrapolated values to know the cancer grade in a neighborhood outsidethe core extracted from the prostate gland during a needle biopsyprocedure. To this end, two models are fitted. The first model is fittedby taking into account the distribution of Gleason grades within theprocessed 2-D or 3-D core. The second model is a pixel-based model,where the pixel value and other features discussed in previousembodiments can be extracted to construct a prediction model. Severalalgorithms, for example the Kalman filtering schemes 26] or similartechniques can be used to make the aforementioned predictions. Once thetwo predictions can be performed, they can be fused to get a uniquecancer grade for the new predicted blocks. An exemplary procedure isshown in FIG. 28.

Referring in more details to the present embodiment, the model of thehistopathology image is given by:

S(m,n)=AS(m−1,n)+Bw(m,n)

Where S is the state vector, A is a state transition matrix, B is aninput matrix, w is random system noise. S represents the noisy observedimage features, which is directly related to the morphology and texturalfeatures extracted for classification purposes. The vector w may bemodeled as a Gaussian noise vector with zero mean, which representvariations in calculated features. As it was mentioned above, pixel- andblock-based model fitting may be performed.

The observation equation is given by:

r(m,n)=h ^(T) S(m,n)+ν(m,n)

Where h is an observation matrix and v is additive noise present inobservations.

The next step is to generate a Kalman filter, for prediction of thecancerous grade of image blocks by estimating the state vector. In orderto execute the 2-D filtering, we divide the 2-D problem into equivalent1-D problems by applying line by line scanning or keeping constant thevariable n or second dimension of the image or block, and thenconstructing a global state vector for all values of n.

For example, the equivalent 1-D Kalman filtering equations for a 1-Dimage scanning process are:

(m)=A

(m−1)

P_(b) (m)=AP_(a)(m−1)A^(T)+BQ_(w)B^(T), Q_(w) is the correlation matrixof noise vector w.

K(m)=P_(b)(m)h^(T) hP_(b)(m)h^(T)+Q_(ν)]⁻¹ Q_(ν) is the correlationmatrix of noise vector v.

(m)=

(m)+K(m)r(m)−h

_(b)(m)]

P _(a)(m)=I−K(m)h]P _(b)(m)

The resulting estimated state vector is used as the input to apreviously trained classifier to obtain the Gleason grade of the 2-Dimage patch. Extensions and generalizations of the theory of estimationby Kalman filtering may be used to generate estimated cancer grade of2-D slices and then reconstruct 3-D images and assign a Gleason gradeusing the classifier for 3-D images developed in previous embodiments.

Cancer Tele-Screening System and Adjudication Schemes (System 9)

A tele-screening system according to an embodiment of the invention isprovided for classification and diagnosis from digitized biopsy imagesin an ideal adjudication scheme. In this embodiment, the final reportwill be produce by integrating pathologists and automatic expert systemsto produce a more accurate cancer diagnosis based on biopsy assessment.FIG. 18 shows the proposed tele-screening system. Several variations ofthe ideal adjudication scheme are addressed in different aspects of thepresent invention.

In one embodiment, an adjudication paradigm for cancer diagnosis biopsyby integrating one expert system is proposed. For instance, theconfiguration of the system is presented in FIG. 19 and includes thefollowing steps: (1) One expert pathologist and one of the computerizedsystems (CAD or expert system) proposed in various embodiments of thisinvention perform independent assessment of all the tissue slidesextracted by a biopsy procedure and produce a diagnostic report andannotated images. (2) An automatic algorithm analyzes the level ofagreement between the two reports including patch analysis of slidescompared to pathologist annotation. The assessment is consideredcomplete if an acceptable level of agreement previously defined isreached. (3) If the minimum acceptable level of agreement is notreached, biopsy samples are assessed by a second pathologist. (4) Theautomatic algorithm of step (2) analyzes the level of agreement amongthe three reports in a patch-based basis and force decision, for exampleby majority voting rule. (5) The misclassified patches are added to thetraining database and the classification model is updated. (6) The finalstep uses cross-validation to ensure that the accuracy and othermeasurement such as sensitivity and specificity remain in an acceptableinterval. For cross-validation, 5-fold, 10-fold, hold-out orleave-one-out methods may be used. This system is able to decrease theinterpretation time of a biopsy and reduce the cost associated to asecond pathologist.

In another embodiment, an adjudication paradigm for cancer diagnosisbiopsy by integrating two independent expert systems is proposed. Theflow diagram of the method is shown in FIG. 20. The configuration of thesystem includes the following steps: (1) One expert pathologist and oneof the computerized systems proposed in various embodiments of thisinvention perform independent assessment of all the tissue slidesextracted by a biopsy procedure and produce a diagnostic report andannotated images. (2) An automatic algorithm analyzes the level ofagreement between the two reports including patch analysis of slidescompared to pathologist annotation. The assessment is consideredcomplete if an acceptable level of agreement previously defined isreached. (3) If the minimum acceptable level of agreement is notreached, biopsy samples are assessed by a second expert system, whichuse different features from the first CAD. For instance, the first CADmay be based on morphological features and the second one on texturalfeatures. (4) The automatic algorithm of step (2) analyzes the level ofagreement among the three reports in a patch-based basis and forcedecision, for example by majority voting rule. (5) The misclassifiedpatches are added to the training database and the classification modelis updated. (6) The final step uses a cross-validation to ensure thatthe accuracy and other measurement such as sensitivity and specificityremain in an acceptable interval. For cross-validation, 5-fold, 10-fold,hold-out or leave-one-out methods may be used. This system is able todecrease the interpretation time of a biopsy and further reduce the costassociated to second and third pathologists.

In another embodiment, an adjudication paradigm for cancer diagnosisbiopsy by integrating three independent expert systems is provided. Theflow diagram of the method is shown in FIG. 21. The configuration of thesystem includes the following steps: (1) Two CAD systems proposed invarious embodiments of this invention perform independent assessment ofall the tissue slides extracted by a biopsy procedure and produce adiagnostic report and annotated images. (2) An automatic algorithmanalyzes the level of agreement between the two reports including patchanalysis of slides compared to pathologist annotation. The assessment isconsidered complete if an acceptable level of agreement previouslydefined is reached. (3) If the minimum acceptable level of agreement isnot reached, biopsy samples are assessed by a second expert system,which use different features from the first CADs. (4) The automaticalgorithm of step (2) analyzes the level of agreement among the threereports in a patch-based basis and force decision, for example bymajority voting rule. This system is able to decrease the interpretationtime of a biopsy from minutes to seconds and further reduce the costassociated to second and third pathologists.

In broad embodiment, the present invention provides several systems for2-D and 3-D classification of histology tissue samples. Such systems useunique feature vectors representing discriminative attributes of eachcancerous pattern, for example Gleason patterns. In addition, othersystems are proposed for cancer prognosis and cancer tele-screening byintegrating one or more of the recognition systems. Moreover, novelalgorithms for image processing such as edge detection, color fractaldimension and color region mapping are provided.

ADDITIONAL EMBODIMENTS

The present invention which is described hereinbefore with reference toflowchart and/or block diagram illustrations of methods, systems,devices, simulations, and computer program products in accordance withsome embodiments of the invention and has been illustrated in detail byusing a computer system. For instance, the flowchart and/or blockdiagrams further illustrate exemplary operations of the computer systemsand methods of FIGS. 1 to 28. It is also conceivable that each block ofthe flowchart and/or block diagram illustrations, and combinations ofblocks in the flowchart and/or block diagram illustrations, may beimplemented by any computer program instructions and/or hardware. Thesecomputer program instructions may be provided to a processor of ageneral purpose computer, a microprocessor, a portable device such ascell phones, a special purpose computer or device, or other programmabledata processing apparatus to produce a device, such that theinstructions, which execute via the processor of the computer or otherprogrammable data processing apparatus, create means for implementingthe functions specified in the flowcharts and/or block diagrams orblocks. The computer program may also be supplied from a remote sourceembodied in a carrier medium such as an electronic signal, including aradio frequency carrier wave or an optical carrier wave.

Further modifications and alternative embodiments of various aspects ofthe invention will be apparent to those skilled in the art in view ofthis description. Accordingly, this description is to be construed asillustrative only and is for the purpose of teaching those skilled inthe art the general manner of carrying out the invention. It is to beunderstood that the forms of the invention shown and described hereinare to be taken as examples of embodiments. Elements and materials maybe substituted for those illustrated and described herein, parts andprocesses may be reversed, and certain features of the invention may beutilized independently, all as would be apparent to one skilled in theart after having the benefit of this description of the invention.Changes may be made in the elements described herein without departingfrom the spirit and scope of the invention as described in the followingclaims.

What is claimed is:
 1. A method, performed by a processing unit, forcomputing pixels along object edges and producing a deinterlaced imagefrom an interlaced source, the method comprising: performing imagefiltering on a collected image depending on the nature of noise in thecollected image; smoothing the filtered image using shape-dependentfilters; calculating gradient vectors in the image using differentkernels; selecting an edge angle; determine threshold values within alocal dynamic range, generating several edge maps, and fusing thegenerated edge maps together.
 2. The method of claim 1, furthercomprising applying an algorithm for image edge-preserving contrastenhancement which is based on HVS, Parameterized Logarithm ImageProcessing operations, and is integrated with morphological log-ratioapproach in order to come up with an effective edge detection operatorthat is sensitive to edges of dark areas in the image.
 3. Methods andsystems for computation of fractal-like dimension for modifying inputdata/image into data describing the image, the method comprising:performing image filtering depending on the nature of noise in theimage; binarizing the image by using the method described in claim 1 or2, or by using an image bit-plane decomposition method; calculating thefractal dimension of binary images by using different grid sizes basedon a Differential Box Counting (DBC) algorithm; and fusing resultingfractal dimensions.
 4. A method for computing the fractal dimension ofcolor images comprising: performing a color model transformation byusing arbitrary functions operating on the original RGB components ofthe image; splitting the image in smaller windows and computing for eachblock the probability of having m points/pixels in a hypercube of thesize of the window; estimating the color fractal dimension by using aweighting function and fitting the curve of logarithm of the windowsizes against the logarithm of the total number of boxes of each windowsize needed to cover the image.
 5. A method for carcinomas color regionmapping comprising: constructing 2-D projections of the RGB color modelof a large variety of biopsy images including normal and canceroustissue; computing 2-D histograms of biopsy images; selecting the mostprominent colors by using an algorithm for local maxima location in 3-Dsurfaces; constructing a color region mapping, in which each colorcorresponds to one class or tissue structure.
 6. A method for imagecolor quantization comprising: performing an image color standardizationprocedure; calculating the distance between each pixel and the colorsconsidered in the color map of claim 4 by using a distance metric;minimizing the calculated distance; and assigning the nearest referencecolor or class to the each pixel.
 7. A textural-based system and methodsfor automatically detecting, classifying, and grading cancerous regionsof a histology image comprising: performing an image colorstandardization procedure; forming texture-based feature vectors byusing spatial and transforms domain image information along with fractalanalysis; selecting the group of features that best describe the images;training a classifier by using the generated feature vectors;classifying histology images according to the Gleason grading system;using the result of classification to determine the Gleason score of theimage; assessing the accuracy of the Gleason grading/scoring system byusing cross-validation methods.
 8. The method of claim 7, whereinforming texture-based feature vectors comprises extraction of a set offeatures from Fourier and wavelet transforms such us wavelet energy andentropy, phase information of coefficients, and spatial domainalgorithms such as statistical spatial filtering, wavelet-based fractaldimension according to the method of claim 3 or others, and gray levelhistogram.
 9. A textural-based system and methods for automaticallydetect, classify, and grade cancerous regions of a histology image(particularly a prostate biopsy images) comprising: performing an imagecolor standardization procedure; describing the image data and formingtexture-based feature vectors by using spatial and transforms domainimage information along with color fractal analysis as claimed in claim4; selecting the group of features that best describe the images;training a classifier by using the generated feature vectors;classifying histology images according to the Gleason grading system;using the result of classification to determine the Gleason score of theimage; assessing the accuracy of the Gleason grading/scoring system byusing cross-validation methods.
 10. The description step of the systemof claim 9, further comprising the extraction of a set of features fromreal, logarithmic, or complex wavelets and multi-wavelets and a newcolor fractal dimension algorithm developed in this invention.Particularly, transform features include joint probability of detailcoefficients, phase information of coefficients, energy distribution ofdetail coefficients, and energy of color channels of the biopsy image. Anew color fractal dimension is also extracted to complete the featurespace.
 11. A morphology- and architectural-based system and methods forautomatically detect, classify, and grade cancerous regions of ahistology image (particularly a prostate biopsy images) comprising:performing an image color standardization procedure; segmenting theimage by applying any proper algorithm including the method proposed inclaims 5 and 6; describing the image data and forming morphology andarchitectural-based feature vectors by extracting the shape, size andarrangement of tissue structures; selecting the group of features thatbest describe the images; training a classifier by using the generatedfeature vectors; classifying cancerous images according to gradingsystems, for example Gleason grading system for prostate carcinomas;using the result of classification to determine the Gleason score of theimage; assessing the accuracy of the Gleason grading/scoring system byusing cross-validation methods.
 12. The description step of the systemof claim 11, further comprising the extraction of features: colorinformation such as color distribution entropy and 2-D color histogramsalong with tissue morphology and architectural features from segmentedtissue structures.
 13. Systems and methods for automatically detecting,classifying and grading cancerous regions of 3-D histological images(particularly a prostate image) comprising: performing an image colorstandardization procedure, segmenting image by applying any properalgorithm including the method proposed in claims 5 and 6; mapping of2-D images or an image into 3-D space by using pre-developed mappingalgorithms, describing the image data and forming color-, texture-,morphology-, and architectural-based feature vectors; training aclassifier within an framework by using generated feature vectors;accessing the cancer grade or scale, for example, according to theGleason grading scale; accessing the most prominent, the second mostprominent pattern and the third prominent pattern in said image data andcomputing a Gleason 3-D score as a sum of Gleason 3-D grades of saidpatterns; accessing the accuracy of Gleason 3-D grading/scoring system;estimating volume of a tumor region.
 14. The systems as claimed in 13,wherein the mapping from 2-D image to 3-D images includes algorithms for3-D reconstruction from a single 2-D image or from several 2-D slides.The processed images nay be gray scaled or colored images.
 15. Thesystems of claims 7, 9, 11 and 13 wherein automatically classifying thebiopsy images is based on machine learning algorithms: lineardiscriminants, Gaussian models, Multivariate Adaptive Regression Splines(MARS), Classification an Regression trees (CART), decision trees,k-nearest neighbor, Bayesian, neural networks, and/or SVM. Boostingalgorithms such as AdaBoost, SVM Boost, etc. can be used to improve theclassification performance.
 16. A method to obtain a multidimensionalGleason grading wherein a revised Gleason value is generated as follows:Revised Gleason Value=2D Gleason Grade+3rd Dimensional Core Grade
 17. Asystem and method for automatically detecting, classifying and gradingcancerous regions from digitized Whole-Slide (particularly a prostateimage) comprising: performing an image color standardization procedure;splitting image into smaller regions; determining and grading a locationof a tumor in the prostate using one or more of the systems of claim 7or 9 or 11 or 13; accessing a cancer stage; accessing the mostprominent, the second most prominent pattern and third prominent patternin said image data and computing a Gleason score as a sum of Gleasongrades of said patterns; generating a cancer map from whole digitizedslides according to Gleason grading/scoring levels associated withslides processing block-elements; accessing the accuracy of Gleasongrading/scoring system; accessing the estimating size of a tumor region.18. The system according to claim 17, wherein said grade estimatingmodule is configured for estimating a size of a tumor in said at leastone tissue region or whole digitized slides.
 19. The system of claim 17,wherein a complete diagnosis is performed for providing usefulinformation for pathologists about the disease severity and thelocalization, size and volume of cancerous tissue.
 20. The system ofclaim 17 wherein complete diagnosing is performed, the outputs of thesystems include: tumor extent in the needle biopsy, number andidentification of positive needle cores, percentage of positive coresfor cancer, length of the tumor in positive cores and length of thecore, percentage of needle core tissue affected by carcinoma and area oftissue positive for cancer, Gleason grade of each selected cancerousregion of interest, percentage of each Gleason pattern in biopsyspecimens, Gleason score of each positive core, including tertiarypattern information, the most frequent Gleason score in slides, numbersof identical Gleason scores in the positive cores, revised Gleason value(claim 16), and localized cancer map including the grade of eachcancerous patch in all needle biopsy slides.
 21. The system of claim 17further comprising methods for surgery quality control, which allow asurgeon to be immediately aware of whether the surrounding tissue of aspecimen is or not positive for cancer, prior the patient is closed up.22. A system for cancer prognosis, wherein the pathology reportincluding all the outputs of systems claimed in 7, 9, 11 and 13 isintegrated with patient information and PSA analysis results toconstruct feature vectors. The resulting feature vectors are used to aidpathologists in cancer prognosis, risk assessment and treatmentplanning.
 23. The system of claim 22, further comprising the use offeature vectors for prostate cancer risk assessment by implementing riskcalculators according to the patient's available information.
 24. Thesystem of claim 22 further comprising a method for cancer predictionbased on 2-D and 3-D histopathology images using cancer map or Gleasongraded block and pixel features to predict the state of the tissue inmissing blocks within the scope of the image or to produce extrapolatedvalues to know the cancer grade in the neighborhoods outside the coreextracted from the prostate gland during a needle biopsy procedure. 25.The method of claim 24 wherein Kalman filtering or similar techniquesare used for estimation and prediction of cancer grade (Gleason grade)outside of the scope of observed histopathology biopsy images.
 26. Asystem and methods for cancer tele-screening and diagnosis based onbiopsy image and an adjudication scheme, wherein expert pathologists andCAD systems are integrated.
 27. The adjudication methods for prostatecancer diagnoses used in the system of claim 26 consist of replacing oneor more pathologists by CAD systems and evaluate the agreement amongthem to produce a faster and less expensive prostate cancer diagnosis.28. The systems of claim 26 wherein said adjudication methods isconfigure for using morphology- or textural-based 2-D classificationsystems for cancer detection and grading, such that the CADs integratedinto the systems perform independent tissue assessment.
 29. Thetele-screening system of claim 26 uses communication networks such asinternet for image and patient data acquisition, pathologists' reportacquisition. The final assessment results are stored in a centralizeddatabase available for pathologist consultation.
 30. A method and systemconstructing an optimal signal detector using a learning algorithm intandem with a prediction algorithm (which may be one in the same) anddata compression algorithm.
 31. A method and system for detectingcarcinomas in a sample image based on compression rates obtained using adetector trained on previously classified carcinoma images as specifiedin claim
 30. 32. A method and system for using detection results tocompute a hard, integer Gleason score of a query image through selectionof a best compressing pattern compressors (PCs) using one of the systemsof claims 30 and
 31. 33. A method and system for using detection resultsto compute a hard, integer Gleason scores of query image pixels andregions through selection of a best compressing pattern compressors(PCs) using one of the systems of claims 30 and
 31. 34. A method andsystem for using detection results to compute a soft, continuous-scaleGleason score of an image (query), comprising: compressing query imagedata using specific pattern PCs and optionally a self PC according tosystems in claims 30 and 31; computing pattern weights determined by: PCcompression rates, PC compression rate statistics, and/or functions ofeach pattern PC compression rates or statistics and query image self PCcompression rates or statistics; combining pattern values and weightsinto a single Gleason score.
 35. A method and system for using detectionresults to compute a soft, continuous-scale Gleason scores of queryimage pixels and regions, comprising: compressing query image pixelsusing specific pattern PCs and optionally a self PC according to systemsin claims 30 and 31; computing pattern weights determined by either bypattern PC compression rate or by comparison of each pattern PCprediction to a query image's self PC prediction on a pixel by pixelbasis; combining pattern values and weights into a single, soft Gleasonscore on a pixel by pixel basis.
 36. A system and method for determiningthe most likely distribution of Gleason patterns within a query image,comprising: the systems of claim in 30, 31 and 33 or 35; determiningpairwise uncertainty statistics between pixel scores or grades; formingthe uncertainty statistics into a matrix; computing the eigenvalues ofsaid matrix; computing the eigenvectors of said matrix; using all orsome of the above features to estimate the most likely currentdistribution of Gleason patterns within the sample; using all or some ofthe above features to estimate the most likely future distribution ofGleason patterns within the sample;
 37. A system and method forautomatically detecting, classifying and grading cancerous regions fromdigitized, Whole-Slides (particularly a prostate image) comprising:segmenting the image into prominent features; optionally downsizing theimage; converting the image into a 1 dimensional signal; compressing thesignal using the systems of claims 30 and 31; assigning Gleason scoresto each pixel using the system of claims 33 or 35, and 36; determining,grading, and scoring a location of a tumor in the prostate using one ormore of the systems of claim 33, 35 or 36; assessing the and grading alocation of a tumor in the prostate using one or more of the systems ofclaim 33, 35 or 36; estimating the current carcinoma patterndistribution of a selected region using the systems of claim 33, 35 or36; estimating the future carcinoma pattern distribution of a selectedregion using the systems of claim 33, 35 or 36; generating a cancer mapfrom whole, digitized slides according to the systems of claim 33, 35 or36; assessing the accuracy of Gleason grading/scoring system;
 38. Thesystem according to claim 37, wherein said grade estimating module isconfigured for estimating a size of a tumor in said at least one tissueregion or whole digitized slides.
 39. The system of claim 37, wherein acomplete diagnosis is performed for providing useful information forpathologists about the disease severity and the localization, size andvolume of cancerous tissue.
 40. The system of claim 37 wherein completediagnosing is performed, the outputs of the systems include: tumorextent in the needle biopsy, number and identification of positiveneedle cores, percentage of positive cores for cancer, length of thetumor in positive cores and length of the core, percentage of needlecore tissue affected by carcinoma and area of tissue positive forcancer, Gleason grade of each selected cancerous region of interest,percentage of each Gleason pattern in biopsy specimens, Gleason score ofeach positive core, including tertiary pattern information, the mostfrequent Gleason score in slides, numbers of identical Gleason scores inthe positive cores, revised Gleason value (claim 16), and localizedcancer map including the grade of each cancerous patch in all needlebiopsy slides.
 41. The system of claim 37 further comprising methods forsurgery quality control, which allow a surgeon to be immediately awareof whether the surrounding tissue of a specimen is or not positive forcancer, prior the patient is closed up.
 42. A system for cancerprognosis, wherein the pathology report including all the outputs ofsystems claimed in 30 through 41 is integrated with patient informationand PSA analysis results to construct feature vectors. The resultingfeature vectors are used to aid pathologists in cancer prognosis, riskassessment and treatment planning.